analysis/TreeMaker/TreeMakerUseFastJet.C

0001 #include <TreeMaker.h>
0002
0003 #include <phool/getClass.h>
0004 #include <phool/PHCompositeNode.h>
0005
0006 // --- calorimeter towers
0007 #include <calobase/RawTower.h>
0008 #include <calobase/RawTowerContainer.h>
0009 #include <calobase/RawTowerGeom.h>
0010 #include <calobase/RawTowerGeomContainer.h>
0011
0012 // --- calorimeter clusters
0013 #include <calobase/RawCluster.h>
0014 #include <calobase/RawClusterv1.h>
0015 #include <calobase/RawClusterContainer.h>
0016
0017 // --- fast jet
0018 #include <fastjet/ClusterSequence.hh>
0019
0020 // --- stuff
0021 #include <TLorentzVector.h>
0022
0023 using std::cout;
0024 using std::endl;
0025 using std::vector;
0026
0027 using fastjet::PseudoJet;
0028 using fastjet::ClusterSequence;
0029 using fastjet::JetDefinition;
0030 using fastjet::antikt_algorithm;
0031
0032
0033
0034 int TreeMaker::UseFastJet(PHCompositeNode *topNode)
0035 {
0036
0037   // -----------------------------------------------------------------------------------------------------
0038   // ---
0039   // -----------------------------------------------------------------------------------------------------
0040
0041   // --- calorimeter tower containers
0042   // RawTowerContainer* cemc_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
0043   // RawTowerContainer* hcalin_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
0044   // RawTowerContainer* hcalout_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
0045
0046   // --- calorimeter geometry objects
0047   // RawTowerGeomContainer* cemc_geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0048   // RawTowerGeomContainer* hcalin_geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0049   // RawTowerGeomContainer* hcalout_geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0050
0051   // --- calorimeter cluster containers
0052   RawClusterContainer* cemc_clusters = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC");
0053   RawClusterContainer* hcalin_clusters = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_HCALIN");
0054   RawClusterContainer* hcalout_clusters = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_HCALOUT");
0055
0056   if ( verbosity > 3 )
0057     {
0058       cout << "cemc_clusters " << cemc_clusters << endl;
0059       cout << "hcalin_clusters " << hcalin_clusters << endl;
0060       cout << "hcalout_clusters " << hcalout_clusters << endl;
0061     }
0062
0063   // --- these nodes are created in CreateNode and filled in CopyAndMakeClusters
0064   RawClusterContainer* new_cemc_clusters = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC_MOD");
0065   RawClusterContainer* new_hcalin_clusters = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_HCALIN_MOD");
0066   RawClusterContainer* new_hcalout_clusters = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_HCALOUT_MOD");
0067
0068   if ( verbosity > 3 )
0069     {
0070       cout << "new_cemc_clusters " << new_cemc_clusters << endl;
0071       cout << "new_hcalin_clusters " << new_hcalin_clusters << endl;
0072       cout << "new_hcalout_clusters " << new_hcalout_clusters << endl;
0073     }
0074
0075   vector<PseudoJet> clusters;
0076   vector<PseudoJet> mod_clusters;
0077   double R = 0.4;
0078   JetDefinition jet_def(antikt_algorithm, R);
0079   // PseudoJet pj_cluster;
0080   // PseudoJet pj_mod_cluster;
0081
0082   // --- loop over cemc clusters
0083   int cluster_counter = 0;
0084   RawClusterContainer::Range cemc_range = cemc_clusters->getClusters();
0085   for ( RawClusterContainer::Iterator cemc_iter = cemc_range.first; cemc_iter != cemc_range.second; ++cemc_iter )
0086     {
0087       // --- get the current cluster
0088       RawCluster* cluster = cemc_iter->second;
0089       double e   = cluster->get_energy();
0090       double phi = cluster->get_phi();
0091       double r   = cluster->get_r();
0092       double z   = cluster->get_z();
0093       double eta = -log(tan(atan2(r,z)/2.0));
0094       double pt  = e/cosh(eta);
0095       TLorentzVector tlv_cluster;
0096       tlv_cluster.SetPtEtaPhiE(pt,eta,phi,e);
0097       clusters.push_back(PseudoJet(tlv_cluster));
0098       ++cluster_counter;
0099     }
0100   ClusterSequence cs_cluster(clusters, jet_def);
0101   vector<PseudoJet> clusterjets = sorted_by_pt(cs_cluster.inclusive_jets());
0102   if ( verbosity > 1 ) cout << "number of clusters " << cluster_counter << " number of cluster jets" << clusterjets.size() << endl;
0103
0104   // --- loop over new_cemc clusters
0105   int mod_cluster_counter = 0;
0106   RawClusterContainer::Range new_cemc_range = new_cemc_clusters->getClusters();
0107   for ( RawClusterContainer::Iterator new_cemc_iter = new_cemc_range.first; new_cemc_iter != new_cemc_range.second; ++new_cemc_iter )
0108     {
0109       // --- get the current cluster
0110       RawCluster* cluster = new_cemc_iter->second;
0111       double e   = cluster->get_energy();
0112       double phi = cluster->get_phi();
0113       double r   = cluster->get_r();
0114       double z   = cluster->get_z();
0115       double eta = -log(tan(atan2(r,z)/2.0));
0116       double pt  = e/cosh(eta);
0117       TLorentzVector tlv_cluster;
0118       tlv_cluster.SetPtEtaPhiE(pt,eta,phi,e);
0119       mod_clusters.push_back(PseudoJet(tlv_cluster));
0120       ++mod_cluster_counter;
0121     }
0122   ClusterSequence cs_mod_cluster(mod_clusters, jet_def);
0123   vector<PseudoJet> mod_clusterjets = sorted_by_pt(cs_mod_cluster.inclusive_jets());
0124   if ( verbosity > 1 ) cout << "number of clusters " << cluster_counter << " number of cluster jets" << mod_clusterjets.size() << endl;
0125
0126
0127   // clusterJet_n = 0;
0128   // for (unsigned i = 0; i < clusterjets.size(); i++) {
0129   //   if (clusterjets[i].pt() < 5) continue;
0130   //   clusterJet_e[i]   = clusterjets[i].e();
0131   //   clusterJet_pt[i]  = clusterjets[i].pt();
0132   //   clusterJet_eta[i] = clusterjets[i].eta();
0133   //   clusterJet_phi[i] = clusterjets[i].phi_std();
0134   //   vector<PseudoJet> constituents = clusterjets[i].constituents();
0135   //   clusterJet_nConst[i] = constituents.size();
0136   //   clusterJet_n++;
0137   // }
0138
0139   // if ( verbosity > 1 )
0140   //   {
0141   //     cout << "process_event: cemc_esum " << cemc_esum << endl;
0142   //     cout << "process_event: new_cemc_esum " << new_cemc_esum << endl;
0143   //   }
0144
0145
0146
0147   return 0;
0148
0149 }