Fork me on GitHub

Changeset 35 in svn for trunk/Utilities


Ignore:
Timestamp:
Nov 17, 2008, 11:07:45 PM (16 years ago)
Author:
Xavier Rouby
Message:

first attemp towers CMS

Location:
trunk/Utilities/Fastjet/plugins/CDFCones
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/Utilities/Fastjet/plugins/CDFCones/interface/CalTower.hh

    r11 r35  
    11#ifndef _CAL_TOWER_HH_
    22#define _CAL_TOWER_HH_
     3
     4// Modified by X. Rouby for integration in Delphes
    35
    46#include <cmath>
    … …  
    810#endif
    911
    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 };
     12const 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
     22const 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}
     27const double TOWER_THETA[ntower] = {
     281.04221324479 ,  1.24151823936 ,  1.4789296564 ,  1.76172836867 ,  2.09858297042 ,  2.49981191976 ,  2.97769357253 ,  3.54683172517 , 
     294.22458468545 ,  5.03156565222 ,  5.99222070913 ,  6.50323511476 ,  8.08260215014 ,  9.38520669791 ,  11.2029938559 ,  13.0015703243 , 
     3014.7731659079 ,  16.5174067026 ,  18.2267592981 ,  19.9103474053 ,  21.6789816306 ,  23.5965023504 ,  25.6731077732 ,  27.9190341722 , 
     3130.3443463769 ,  32.9586699802 ,  35.7708584993 ,  38.7885912317 ,  42.017901925 ,  45.4626450521 ,  49.1239157061 ,  52.9994507062 , 
     3257.0830516071 ,  61.3640831827 ,  65.8271109013 ,  70.4517445074 ,  75.2127486134 ,  80.0804628398 ,  85.02154355, 90
     33};
    1234
    1335class CalTower
    1436{
    1537 public:
    16 
    1738  double Et,eta,phi;
    1839  int iEta,iPhi;
    1940
    2041  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)
    2243  {
    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
    2748            iEta = 4 + i;
    2849            break;
    2950          }
    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
    3455            iEta = 47 - i;
    3556            break;
    3657          }
    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;
    4061      else
    41         iPhi = int(phi/2/M_PI*24)%24;
     62        iPhi = int(phi/(2*pi)*24)%24;
    4263    }
    43     else{
     64    else{ // particle beyond detector reach :
    4465      iEta = -1;
    4566      iPhi = -1;
    4667    }
    4768  }
    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) {}
    4970  CalTower(const CalTower& c): Et(c.Et), eta(c.eta), phi(c.phi), iEta(c.iEta), iPhi(c.iPhi) {}
    50   bool isEqual(CalTower c)
     71  bool isEqual(const CalTower& c)
    5172  {
    5273    return Et == c.Et && eta == c.eta && phi == c.phi && iEta == c.iEta && iPhi == c.iPhi;
    5374  }
     75
     76// converts theta angle [degrees] into pseudorapidity
     77 double pseudorapidity(const double theta) {
     78         return -log(tan( (theta/2.) * (pi/180.) ));
     79         }
     80
    5481};
    5582
     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
     10110   30.0     -1.19999982651       14
     10211   33.524   -1.10000635195       15
     10312   36.822   -1.00361552815       16
     10413   40.261   -0.916007799535      17
     10514   43.614   -0.822472442947      18
     10615   47.436   -0.72264587494       19
     10716   51.79    -0.616250691646      20
     10817   56.735   -0.503273260393      21
     10918   62.31    -0.384075299436      22
     11019   68.516   -0.259479460739      23
     11120   75.297   -0.130817437415      24
     11221   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
     12710   30.0     -1.19999982651      37
     12811   33.524   -1.10000635195      36
     12912   36.822   -1.00361552815      35
     13013   40.261   -0.916007799535     34
     13114   43.614   -0.822472442947     33
     13215   47.436   -0.72264587494      32
     13316   51.79    -0.616250691646     31
     13417   56.735   -0.503273260393     30
     13518   62.31    -0.384075299436     29
     13619   68.516   -0.259479460739     28
     13720   75.297   -0.130817437415     27
     13821   82.526   1.11022302463e-16   26
     139*****/
     140
    56141#endif
     142

  • trunk/Utilities/Fastjet/plugins/CDFCones/interface/Centroid.hh

    r11 r35  
    1111{
    1212 public:
    13 
    1413  double Et,eta,phi;
    1514
    … …  
    1716  Centroid(double centroidEt, double centroidEta, double centroidPhi): Et(centroidEt), eta(centroidEta), phi(centroidPhi) {}
    1817  Centroid(const Centroid& c): Et(c.Et), eta(c.eta), phi(c.phi) {}
    19   void add(Centroid c)
     18  void add(const Centroid& c)
    2019  {
    2120    double newEt = Et + c.Et;
    … …  
    3332    Et = newEt;
    3433  }
    35   void subtract(Centroid c)
     34  void subtract(const Centroid& c)
    3635  {
    3736    double newEt = Et - c.Et;
    … …  
    4948    Et = newEt;
    5049  }
    51   bool isEqual(Centroid c)
     50  bool isEqual(const Centroid& c)
    5251  {
    5352    return Et == c.Et && eta == c.eta && phi == c.phi;

  • trunk/Utilities/Fastjet/plugins/CDFCones/interface/Cluster.hh

    r11 r35  
    3030    pt_tilde = 0.0;
    3131  }
    32   void addTower(PhysicsTower p)
     32  void addTower(const PhysicsTower& p)
    3333  {
    3434    towerList.push_back(p);
    … …  
    3737    pt_tilde += p.fourVector.pt();
    3838  }
    39   void removeTower(PhysicsTower p)
     39  void removeTower(const PhysicsTower& p)
    4040  {
    4141    for(std::vector<PhysicsTower>::iterator towerIter = towerList.begin(); towerIter != towerList.end(); towerIter++)

  • trunk/Utilities/Fastjet/plugins/CDFCones/interface/LorentzVector.hh

    r11 r35  
    3030    return r;
    3131  }
    32   void add(LorentzVector v)
     32  void add(const LorentzVector& v)
    3333  {
    3434    px += v.px;
    … …  
    3737    E  += v.E;
    3838  }
    39   void subtract(LorentzVector v)
     39  void subtract(const LorentzVector& v)
    4040  {
    4141    px -= v.px;
    … …  
    4444    E  -= v.E;
    4545  }
    46   bool isEqual(LorentzVector v)
     46  bool isEqual(const LorentzVector& v)
    4747  {
    4848    return px == v.px && py == v.py && pz == v.pz && E == v.E;

  • trunk/Utilities/Fastjet/plugins/CDFCones/interface/MidPointAlgorithm.hh

    r11 r35  
    1515  double _coneRadius;
    1616  double _coneAreaFraction;
    17   int    _maxPairSize;
     17  unsigned int    _maxPairSize;
    1818  int    _maxIterations;
    1919  double _overlapThreshold;
    … …  
    4444                   std::vector<Cluster>& stableCones, bool reduceConeSize);
    4545  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);
    4747  void splitAndMerge(std::vector<Cluster>& stableCones, std::vector<Cluster>& jets);
    4848  void run(std::vector<PhysicsTower>& particles, std::vector<Cluster>& jets);

  • trunk/Utilities/Fastjet/plugins/CDFCones/interface/PhysicsTower.hh

    r11 r35  
    88{
    99 public:
    10 
    1110  LorentzVector fourVector;
    1211  CalTower calTower;
     12  /// addition by GPS (2008-08-15) for tracking within fastjet
     13  int fjindex;
    1314
    1415  PhysicsTower(): fourVector(LorentzVector()), calTower(CalTower()), fjindex(-1) {}
    15   PhysicsTower(LorentzVector v, CalTower c): fourVector(v), calTower(c), fjindex(-1) {}
     16  PhysicsTower(const LorentzVector& v, const CalTower& c): fourVector(v), calTower(c), fjindex(-1) {}
    1617  PhysicsTower(const PhysicsTower& p): fourVector(p.fourVector), calTower(p.calTower), fjindex(p.fjindex) {}
    17   PhysicsTower(CalTower c):
     18  PhysicsTower(const CalTower& c):
    1819    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(LorentzVector v): 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) {}
    2021  double Et()   const {return calTower.Et;}
    2122  double eta()  const {return calTower.eta;}
    … …  
    2324  int    iEta() const {return calTower.iEta;}
    2425  int    iPhi() const {return calTower.iPhi;}
    25   bool isEqual(PhysicsTower p)
     26  bool isEqual(const PhysicsTower& p)
    2627  {
    2728    return fourVector.isEqual(p.fourVector) && calTower.isEqual(p.calTower);
    2829  }
    29   /// addition by GPS (2008-08-15) for tracking within fastjet
    30   int fjindex;
     30
    3131};
    3232

  • trunk/Utilities/Fastjet/plugins/CDFCones/src/MidPointAlgorithm.cc

    r11 r35  
    1818  std::vector< std::vector<bool> > distanceOK;
    1919  distanceOK.resize(stableCones.size() - 1);
    20   for(int nCluster1 = 1; nCluster1 < stableCones.size(); nCluster1++){
     20  for(unsigned int nCluster1 = 1; nCluster1 < stableCones.size(); nCluster1++){
    2121    distanceOK[nCluster1 - 1].resize(nCluster1);
    2222    double cluster1Rapidity = stableCones[nCluster1].fourVector.y();
    2323    double cluster1Phi      = stableCones[nCluster1].fourVector.phi();
    24     for(int nCluster2 = 0; nCluster2 < nCluster1; nCluster2++){
     24    for(unsigned int nCluster2 = 0; nCluster2 < nCluster1; nCluster2++){
    2525      double cluster2Rapidity = stableCones[nCluster2].fourVector.y();
    2626      double cluster2Phi      = stableCones[nCluster2].fourVector.phi();
    … …  
    3737  std::vector< std::vector<int> > pairs(0);
    3838  std::vector<int> testPair(0);
    39   int maxClustersInPair = _maxPairSize;
     39  unsigned int maxClustersInPair = _maxPairSize;
    4040  if(!maxClustersInPair)
    4141    maxClustersInPair = stableCones.size();
    … …  
    4444  // Loop over all combinations. Calculate MidPoint. Make midPointClusters.
    4545  bool reduceConeSize = false;
    46   for(int iPair = 0; iPair < pairs.size(); iPair++){
     46  for(unsigned int iPair = 0; iPair < pairs.size(); iPair++){
    4747    // Calculate rapidity, phi and pT of MidPoint.
    4848    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++)
    5050      midPoint.add(stableCones[pairs[iPair][iPairMember]].fourVector);
    5151    iterateCone(midPoint.y(),midPoint.phi(),midPoint.pt(),towers,stableCones,reduceConeSize);
    … …  
    116116
    117117void 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)
    119119{
    120120  // Recursively adds clusters to pairs, triplets, ... whose mid-points are then calculated.
    … …  
    124124  if(testPair.size())
    125125    nextClusterStart = testPair.back() + 1;
    126   for(int nextCluster = nextClusterStart; nextCluster <= distanceOK.size(); nextCluster++){
     126  for(unsigned int nextCluster = nextClusterStart; nextCluster <= distanceOK.size(); nextCluster++){
    127127    // Is new SeedCone less than 2*_coneRadius apart from all clusters in testPair?
    128128    bool addCluster = true;
    129     for(int iCluster = 0; iCluster < testPair.size() && addCluster; iCluster++)
     129    for(unsigned int iCluster = 0; iCluster < testPair.size() && addCluster; iCluster++)
    130130      if(!distanceOK[nextCluster - 1][testPair[iCluster]])
    131131        addCluster = false;
Note: See TracChangeset for help on using the changeset viewer.