analysis/SimpleTrackingAnalysis/STACalorimeterCharacterization.C

0001 #include "STACalorimeterCharacterization.h"
0002
0003 #include <phool/getClass.h>
0004 #include <fun4all/Fun4AllServer.h>
0005 #include <fun4all/Fun4AllReturnCodes.h>
0006
0007 #include <phool/PHCompositeNode.h>
0008 //#include <g4hough/PHG4HoughTransform.h> // CAUSES ERRORS DUE TO VertexFinder.h FAILING TO FIND <Eigen/LU>
0009
0010
0011 #include <g4main/PHG4TruthInfoContainer.h>
0012 #include <g4main/PHG4Particle.h>
0013 #include <g4main/PHG4VtxPoint.h>
0014
0015 #include <g4hough/SvtxVertexMap.h>
0016 #include <g4hough/SvtxVertex.h>
0017 #include <g4hough/SvtxTrackMap.h>
0018 #include <g4hough/SvtxTrack.h>
0019
0020 #include <g4eval/SvtxEvalStack.h>
0021 #include <g4eval/SvtxTrackEval.h>
0022 #include <g4eval/SvtxVertexEval.h>
0023 #include <g4eval/SvtxTruthEval.h>
0024
0025 // --- common to all calorimeters
0026 #include <calobase/RawTowerGeomContainer.h>
0027 #include <calobase/RawTowerContainer.h>
0028 #include <calobase/RawTower.h>
0029 #include <calobase/RawCluster.h>
0030 #include <calobase/RawClusterContainer.h>
0031 #include <g4eval/CaloEvalStack.h>
0032 #include <g4eval/CaloRawClusterEval.h>
0033 #include <g4eval/CaloRawTowerEval.h>
0034
0035 #include <TH1D.h>
0036 #include <TH2D.h>
0037 #include <TProfile2D.h>
0038
0039 #include <iostream>
0040 //#include <cassert>
0041
0042 using namespace std;
0043
0044
0045
0046 STACalorimeterCharacterization::STACalorimeterCharacterization(const string &name) : SubsysReco(name)
0047 {
0048   cout << "Class constructor called " << endl;
0049   nevents = 0;
0050   nerrors = 0;
0051   nwarnings = 0;
0052   nlayers = 7;
0053   verbosity = 0;
0054   magneticfield = 1.4; // default is 1.5 but Mike's tracking stuff is using 1.4 for now...
0055   //docalocuts = false;
0056 }
0057
0058
0059
0060 int STACalorimeterCharacterization::Init(PHCompositeNode *topNode)
0061 {
0062
0063   // --------------------------------------------------------------------------------
0064   // --- This is the class initialization method
0065   // --- Get a pointer to the Fun4All server to register histograms
0066   // --- Histograms are declared in class header file
0067   // --- Histograms are instantiated here and then registered with the Fun4All server
0068   // --------------------------------------------------------------------------------
0069
0070
0071   // --- get instance of
0072   Fun4AllServer *se = Fun4AllServer::instance();
0073
0074
0075
0076
0077
0078
0079   // --- some basic calorimeter performance histograms
0080
0081   _energy_ratio_emc = new TH2D("energy_ratio_emc", "", 300,0.0,30.0, 100,0.0,2.0);
0082   _energy_ratio_hci = new TH2D("energy_ratio_hci", "", 300,0.0,30.0, 100,0.0,2.0);
0083   _energy_ratio_hco = new TH2D("energy_ratio_hco", "", 300,0.0,30.0, 100,0.0,2.0);
0084   _energy_ratio_hct = new TH2D("energy_ratio_hct", "", 300,0.0,30.0, 100,0.0,2.0);
0085   _energy_ratio_tot_dumb = new TH2D("energy_ratio_tot_dumb", "", 300,0.0,30.0, 100,0.0,2.0);
0086   _energy_ratio_tot_smart = new TH2D("energy_ratio_tot_smart", "", 300,0.0,30.0, 100,0.0,2.0);
0087   se->registerHisto(_energy_ratio_emc);
0088   se->registerHisto(_energy_ratio_hci);
0089   se->registerHisto(_energy_ratio_hco);
0090   se->registerHisto(_energy_ratio_hct);
0091   se->registerHisto(_energy_ratio_tot_dumb);
0092   se->registerHisto(_energy_ratio_tot_smart);
0093
0094   _energy_ratio_elb_emc = new TH2D("energy_ratio_elb_emc", "", 300,0.0,30.0, 100,0.0,2.0);
0095   _energy_ratio_elb_hci = new TH2D("energy_ratio_elb_hci", "", 300,0.0,30.0, 100,0.0,2.0);
0096   _energy_ratio_elb_hco = new TH2D("energy_ratio_elb_hco", "", 300,0.0,30.0, 100,0.0,2.0);
0097   _energy_ratio_elb_hct = new TH2D("energy_ratio_elb_hct", "", 300,0.0,30.0, 100,0.0,2.0);
0098   _energy_ratio_elb_tot_dumb = new TH2D("energy_ratio_elb_tot_dumb", "", 300,0.0,30.0, 100,0.0,2.0);
0099   _energy_ratio_elb_tot_smart = new TH2D("energy_ratio_elb_tot_smart", "", 300,0.0,30.0, 100,0.0,2.0);
0100   se->registerHisto(_energy_ratio_elb_emc);
0101   se->registerHisto(_energy_ratio_elb_hci);
0102   se->registerHisto(_energy_ratio_elb_hco);
0103   se->registerHisto(_energy_ratio_elb_hct);
0104   se->registerHisto(_energy_ratio_elb_tot_dumb);
0105   se->registerHisto(_energy_ratio_elb_tot_smart);
0106
0107   _energy_ratio_eub_emc = new TH2D("energy_ratio_eub_emc", "", 300,0.0,30.0, 100,0.0,2.0);
0108   _energy_ratio_eub_hci = new TH2D("energy_ratio_eub_hci", "", 300,0.0,30.0, 100,0.0,2.0);
0109   _energy_ratio_eub_hco = new TH2D("energy_ratio_eub_hco", "", 300,0.0,30.0, 100,0.0,2.0);
0110   _energy_ratio_eub_hct = new TH2D("energy_ratio_eub_hct", "", 300,0.0,30.0, 100,0.0,2.0);
0111   _energy_ratio_eub_tot_dumb = new TH2D("energy_ratio_eub_tot_dumb", "", 300,0.0,30.0, 100,0.0,2.0);
0112   _energy_ratio_eub_tot_smart = new TH2D("energy_ratio_eub_tot_smart", "", 300,0.0,30.0, 100,0.0,2.0);
0113   se->registerHisto(_energy_ratio_eub_emc);
0114   se->registerHisto(_energy_ratio_eub_hci);
0115   se->registerHisto(_energy_ratio_eub_hco);
0116   se->registerHisto(_energy_ratio_eub_hct);
0117   se->registerHisto(_energy_ratio_eub_tot_dumb);
0118   se->registerHisto(_energy_ratio_eub_tot_smart);
0119
0120   _energy_dphi_emc = new TH2D("energy_dphi_emc", "", 300,0.0,30.0, 100,-1.0,1.0);
0121   _energy_dphi_hci = new TH2D("energy_dphi_hci", "", 300,0.0,30.0, 100,-1.0,1.0);
0122   _energy_dphi_hco = new TH2D("energy_dphi_hco", "", 300,0.0,30.0, 100,-1.0,1.0);
0123   se->registerHisto(_energy_dphi_emc);
0124   se->registerHisto(_energy_dphi_hci);
0125   se->registerHisto(_energy_dphi_hco);
0126
0127   _energy_deta_emc = new TH2D("energy_deta_emc", "", 300,0.0,30.0, 100,-1.0,1.0);
0128   _energy_deta_hci = new TH2D("energy_deta_hci", "", 300,0.0,30.0, 100,-1.0,1.0);
0129   _energy_deta_hco = new TH2D("energy_deta_hco", "", 300,0.0,30.0, 100,-1.0,1.0);
0130   se->registerHisto(_energy_deta_emc);
0131   se->registerHisto(_energy_deta_hci);
0132   se->registerHisto(_energy_deta_hco);
0133
0134
0135
0136   _towers_3x3_emc = new TProfile2D("towers_3x3_emc", "", 3,-1.5,1.5, 3,-1.5,1.5, 0.0,9998.0,"");
0137   _towers_5x5_emc = new TProfile2D("towers_5x5_emc", "", 5,-2.5,2.5, 5,-2.5,2.5, 0.0,9998.0,"");
0138   _towers_7x7_emc = new TProfile2D("towers_7x7_emc", "", 7,-3.5,3.5, 7,-3.5,3.5, 0.0,9998.0,"");
0139   _towers_9x9_emc = new TProfile2D("towers_9x9_emc", "", 9,-4.5,4.5, 9,-4.5,4.5, 0.0,9998.0,"");
0140   se->registerHisto(_towers_3x3_emc);
0141   se->registerHisto(_towers_5x5_emc);
0142   se->registerHisto(_towers_7x7_emc);
0143   se->registerHisto(_towers_9x9_emc);
0144
0145   for ( int i = 0; i < 10; ++i )
0146     {
0147       _towersum_energy_emc[i] = new TH2D(Form("towersum_energy_emc_%d",i), "", 300,0.0,30.0, 100,0.0,2.0);
0148       _tower_energy_emc[i] = new TH2D(Form("tower_energy_emc_%d",i), "", 300,0.0,30.0, 100,0.0,2.0);
0149       se->registerHisto(_towersum_energy_emc[i]);
0150       se->registerHisto(_tower_energy_emc[i]);
0151     }
0152
0153   _towers_3x3_hci = new TProfile2D("towers_3x3_hci", "", 3,-1.5,1.5, 3,-1.5,1.5, 0.0,9998.0,"");
0154   _towers_5x5_hci = new TProfile2D("towers_5x5_hci", "", 5,-2.5,2.5, 5,-2.5,2.5, 0.0,9998.0,"");
0155   _towers_7x7_hci = new TProfile2D("towers_7x7_hci", "", 7,-3.5,3.5, 7,-3.5,3.5, 0.0,9998.0,"");
0156   _towers_9x9_hci = new TProfile2D("towers_9x9_hci", "", 9,-4.5,4.5, 9,-4.5,4.5, 0.0,9998.0,"");
0157   se->registerHisto(_towers_3x3_hci);
0158   se->registerHisto(_towers_5x5_hci);
0159   se->registerHisto(_towers_7x7_hci);
0160   se->registerHisto(_towers_9x9_hci);
0161
0162   for ( int i = 0; i < 10; ++i )
0163     {
0164       _towersum_energy_hci[i] = new TH2D(Form("towersum_energy_hci_%d",i), "", 300,0.0,30.0, 100,0.0,2.0);
0165       _tower_energy_hci[i] = new TH2D(Form("tower_energy_hci_%d",i), "", 300,0.0,30.0, 100,0.0,2.0);
0166       se->registerHisto(_towersum_energy_hci[i]);
0167       se->registerHisto(_tower_energy_hci[i]);
0168     }
0169
0170   _towers_3x3_hco = new TProfile2D("towers_3x3_hco", "", 3,-1.5,1.5, 3,-1.5,1.5, 0.0,9998.0,"");
0171   _towers_5x5_hco = new TProfile2D("towers_5x5_hco", "", 5,-2.5,2.5, 5,-2.5,2.5, 0.0,9998.0,"");
0172   _towers_7x7_hco = new TProfile2D("towers_7x7_hco", "", 7,-3.5,3.5, 7,-3.5,3.5, 0.0,9998.0,"");
0173   _towers_9x9_hco = new TProfile2D("towers_9x9_hco", "", 9,-4.5,4.5, 9,-4.5,4.5, 0.0,9998.0,"");
0174
0175   for ( int i = 0; i < 10; ++i )
0176     {
0177       _towersum_energy_hco[i] = new TH2D(Form("towersum_energy_hco_%d",i), "", 300,0.0,30.0, 100,0.0,2.0);
0178       _tower_energy_hco[i] = new TH2D(Form("tower_energy_hco_%d",i), "", 300,0.0,30.0, 100,0.0,2.0);
0179       se->registerHisto(_towersum_energy_hco[i]);
0180       se->registerHisto(_tower_energy_hco[i]);
0181     }
0182   se->registerHisto(_towers_3x3_hco);
0183   se->registerHisto(_towers_5x5_hco);
0184   se->registerHisto(_towers_7x7_hco);
0185   se->registerHisto(_towers_9x9_hco);
0186
0187
0188
0189   return 0;
0190
0191 }
0192
0193
0194
0195 int STACalorimeterCharacterization::process_event(PHCompositeNode *topNode)
0196 {
0197
0198   // --- This is the class process_event method
0199   // --- This is where the bulk of the analysis is done
0200   // --- Here we get the various data nodes we need to do the analysis
0201   // --- Then we use variables (accessed through class methods) to perform calculations
0202
0203   if ( verbosity > -1 )
0204     {
0205       cout << endl;
0206       cout << "------------------------------------------------------------------------------------" << endl;
0207       cout << "Now processing event number " << nevents << endl; // would be good to add verbosity switch
0208     }
0209
0210   ++nevents; // You may as youtself, why ++nevents (pre-increment) rather
0211   // than nevents++ (post-increment)?  The short answer is performance.
0212   // For simple types it probably doesn't matter, but it can really help
0213   // for complex types (like the iterators below).
0214
0215
0216   // --- Truth level information
0217   PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,"G4TruthInfo");
0218   if ( !truthinfo )
0219     {
0220       cerr << PHWHERE << " ERROR: Can't find G4TruthInfo" << endl;
0221       exit(-1);
0222     }
0223
0224   // --- SvtxTrackMap
0225   SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
0226   if ( !trackmap )
0227     {
0228       cerr << PHWHERE << " ERROR: Can't find SvtxTrackMap" << endl;
0229       exit(-1);
0230     }
0231
0232   // --- SvtxVertexMap
0233   SvtxVertexMap* vertexmap = findNode::getClass<SvtxVertexMap>(topNode,"SvtxVertexMap");
0234   if ( !vertexmap )
0235     {
0236       cerr << PHWHERE << " ERROR: Can't find SvtxVertexMap" << endl;
0237       exit(-1);
0238     }
0239
0240   // --- Raw cluster classes
0241   bool clusters_available = true;
0242   RawClusterContainer *emc_clustercontainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC");
0243   RawClusterContainer *hci_clustercontainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_HCALIN");
0244   RawClusterContainer *hco_clustercontainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_HCALOUT");
0245   if ( !emc_clustercontainer || !hci_clustercontainer || !hco_clustercontainer )
0246     {
0247       if ( verbosity > -1 )
0248     {
0249       cerr << PHWHERE << " WARNING: Can't find cluster nodes" << endl;
0250       cerr << PHWHERE << "  emc_clustercontainer " << emc_clustercontainer << endl;
0251       cerr << PHWHERE << "  hci_clustercontainer " << hci_clustercontainer << endl;
0252       cerr << PHWHERE << "  hco_clustercontainer " << hco_clustercontainer << endl;
0253     }
0254       clusters_available = false;
0255       ++nwarnings;
0256     }
0257
0258   // --- Tower geometry
0259   RawTowerGeomContainer *emc_towergeo = findNode::getClass<RawTowerGeomContainer>(topNode,"TOWERGEOM_CEMC");
0260   RawTowerGeomContainer *hci_towergeo = findNode::getClass<RawTowerGeomContainer>(topNode,"TOWERGEOM_HCALIN");
0261   RawTowerGeomContainer *hco_towergeo = findNode::getClass<RawTowerGeomContainer>(topNode,"TOWERGEOM_HCALOUT");
0262   if ( !emc_towergeo || !hci_towergeo || !hco_towergeo )
0263     {
0264       if ( verbosity > -1 )
0265     {
0266       cerr << PHWHERE << " WARNING: Can't find tower geometry nodes" << endl;
0267       cerr << PHWHERE << "  emc_towergeo " << emc_towergeo << endl;
0268       cerr << PHWHERE << "  hci_towergeo " << hci_towergeo << endl;
0269       cerr << PHWHERE << "  hco_towergeo " << hco_towergeo << endl;
0270     }
0271       ++nwarnings;
0272     }
0273
0274   // --- Tower container
0275   RawTowerContainer *emc_towercontainer = findNode::getClass<RawTowerContainer>(topNode,"TOWER_CALIB_CEMC");
0276   RawTowerContainer *hci_towercontainer = findNode::getClass<RawTowerContainer>(topNode,"TOWER_CALIB_HCALIN");
0277   RawTowerContainer *hco_towercontainer = findNode::getClass<RawTowerContainer>(topNode,"TOWER_CALIB_HCALOUT");
0278   if ( !emc_towercontainer || !hci_towercontainer || !hco_towercontainer )
0279     {
0280       if ( verbosity > -1 )
0281     {
0282       cerr << PHWHERE << " WARNING: Can't find tower container nodes" << endl;
0283       cerr << PHWHERE << "  emc_towercontainer " << emc_towercontainer << endl;
0284       cerr << PHWHERE << "  hci_towercontainer " << hci_towercontainer << endl;
0285       cerr << PHWHERE << "  hco_towercontainer " << hco_towercontainer << endl;
0286     }
0287       ++nwarnings;
0288     }
0289
0290
0291
0292
0293   // --- step 1: loop over all possible towers and make a map of energy, tower address
0294   // --- step 2: loop over map in reverse order and fill a vector of tower addresses
0295   // --- step 3: use vector of tower addresses (which are ordered by energy highest to lowest) as desired
0296   // --- note: the vector is not really needed, as one can just get anything that's needed from the map
0297   // ---       itself, but personally I like vectors
0298
0299   //int nphi = emc_towergeo->get_phibins();
0300   //int neta = emc_towergeo->get_etabins();
0301
0302
0303
0304   vector<RawCluster*> emc_clusters = get_ordered_clusters(emc_clustercontainer);
0305   vector<RawCluster*> hci_clusters = get_ordered_clusters(hci_clustercontainer);
0306   vector<RawCluster*> hco_clusters = get_ordered_clusters(hco_clustercontainer);
0307
0308
0309
0310   // --- Create SVTX eval stack
0311   SvtxEvalStack svtxevalstack(topNode);
0312   // --- Get evaluator objects from the eval stack
0313   //SvtxVertexEval*   vertexeval = svtxevalstack.get_vertex_eval();
0314   SvtxTrackEval*     trackeval = svtxevalstack.get_track_eval();
0315   SvtxTruthEval*     trutheval = svtxevalstack.get_truth_eval();
0316
0317   // --- Create calorimter eval stacks
0318   CaloEvalStack emc_caloevalstack(topNode,"CLUSTER_CEMC");
0319   CaloEvalStack hci_caloevalstack(topNode,"CLUSTER_HCALIN");
0320   CaloEvalStack hco_caloevalstack(topNode,"CLUSTER_HCALOUT");
0321   // --- get the evaluator objects
0322   CaloRawClusterEval *emc_rawclustereval = emc_caloevalstack.get_rawcluster_eval();
0323   CaloRawClusterEval *hci_rawclustereval = hci_caloevalstack.get_rawcluster_eval();
0324   CaloRawClusterEval *hco_rawclustereval = hco_caloevalstack.get_rawcluster_eval();
0325   emc_rawclustereval->set_verbosity(0); // temp while resolving issues
0326   hci_rawclustereval->set_verbosity(0); // temp while resolving issues
0327   hco_rawclustereval->set_verbosity(0); // temp while resolving issues
0328
0329   if ( verbosity > 0 || ( !clusters_available && verbosity > -1 ) )
0330     {
0331       cout << "Eval stack memory addresses..." << endl;
0332       cout << &svtxevalstack << endl;
0333       cout << &emc_caloevalstack << endl;
0334       cout << &hci_caloevalstack << endl;
0335       cout << &hco_caloevalstack << endl;
0336       cout << "Evaluator addresses..." << endl;
0337       cout << trackeval << endl;
0338       cout << emc_rawclustereval << endl;
0339       cout << hci_rawclustereval << endl;
0340       cout << hco_rawclustereval << endl;
0341     }
0342
0343   if ( verbosity > 0 ) cout << "Now going to loop over truth partcles..." << endl; // need verbosity switch
0344
0345   // --- Loop over all truth particles
0346   PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
0347   for ( PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter )
0348     {
0349       PHG4Particle* g4particle = iter->second; // You may ask yourself, why second?
0350       // In C++ the iterator is a map, which has two members
0351       // first is the key (analogous the index of an arry),
0352       // second is the value (analogous to the value stored for the array index)
0353       int particleID = g4particle->get_pid();
0354
0355       if ( trutheval->get_embed(g4particle) <= 0 && fabs(particleID) == 11 && verbosity > 0 )
0356         {
0357           cout << "NON EMBEDDED ELECTRON!!!  WHEE!!! " << particleID << " " << iter->first << endl;
0358         }
0359
0360       if ( trutheval->get_embed(g4particle) <= 0 ) continue; // only look at embedded particles // no good for hits files
0361       bool iselectron = fabs(particleID) == 11;
0362       bool ispion = fabs(particleID) == 211;
0363       if ( verbosity > 0 ) cout << "embedded particle ID is " << particleID << " ispion " << ispion << " iselectron " << iselectron << " " << iter->first << endl;
0364
0365       set<PHG4Hit*> g4hits = trutheval->all_truth_hits(g4particle);
0366       //float ng4hits = g4hits.size();
0367
0368       float truept = sqrt(pow(g4particle->get_px(),2)+pow(g4particle->get_py(),2));
0369       float true_energy = g4particle->get_e();
0370
0371       // --- Get the reconsructed SvtxTrack based on the best candidate from the truth info
0372       SvtxTrack* track = trackeval->best_track_from(g4particle);
0373       if (!track) continue;
0374       float recopt = track->get_pt();
0375       //float recop = track->get_p();
0376
0377       if ( verbosity > 0 )
0378     {
0379       cout << "truept is " << truept << endl;
0380       cout << "recopt is " << recopt << endl;
0381       cout << "true energy is " << true_energy << endl;
0382     }
0383
0384       // ----------------------------------------------------------------------
0385       // ----------------------------------------------------------------------
0386       // ----------------------------------------------------------------------
0387
0388       // --- Get the clusters from the best candidate from the truth info using the eval
0389       RawCluster* emc_bestcluster = NULL;
0390       RawCluster* hci_bestcluster = NULL;
0391       RawCluster* hco_bestcluster = NULL;
0392       if ( emc_rawclustereval ) emc_bestcluster = emc_rawclustereval->best_cluster_from(g4particle);
0393       if ( hci_rawclustereval ) hci_bestcluster = hci_rawclustereval->best_cluster_from(g4particle);
0394       if ( hco_rawclustereval ) hco_bestcluster = hco_rawclustereval->best_cluster_from(g4particle);
0395
0396       if ( verbosity > 5 )
0397     {
0398       cout << "Cluster addresses from best cluster " << endl;
0399       cout << emc_bestcluster << endl;
0400       cout << hci_bestcluster << endl;
0401       cout << hco_bestcluster << endl;
0402     }
0403
0404       // --- If that didn't work, take the largest cluster in the event
0405       // --- This is terrible for more than one particle, so I need to
0406       // --- develop my own track-cluster association...
0407       if ( !emc_bestcluster && emc_clusters.size() ) emc_bestcluster = emc_clusters[0];
0408       if ( !hci_bestcluster && hci_clusters.size() ) hci_bestcluster = hci_clusters[0];
0409       if ( !hco_bestcluster && hco_clusters.size() ) hco_bestcluster = hco_clusters[0];
0410
0411       if ( verbosity > 5 )
0412     {
0413       cout << "Cluster addresses from my (very bad) association " << endl;
0414       cout << emc_bestcluster << endl;
0415       cout << hci_bestcluster << endl;
0416       cout << hco_bestcluster << endl;
0417     }
0418
0419       // --- Get a vector of the towers from the cluster
0420       vector<RawTower*> emc_towers_from_bestcluster = get_ordered_towers(emc_bestcluster,emc_towercontainer);
0421       vector<RawTower*> hci_towers_from_bestcluster = get_ordered_towers(hci_bestcluster,hci_towercontainer);
0422       vector<RawTower*> hco_towers_from_bestcluster = get_ordered_towers(hco_bestcluster,hco_towercontainer);
0423
0424       // --- Inspect the towers (fills a bunch of histograms, prints to screen if verbose)
0425       inspect_ordered_towers(emc_towers_from_bestcluster,true_energy,SvtxTrack::CEMC);
0426       inspect_ordered_towers(hci_towers_from_bestcluster,true_energy,SvtxTrack::HCALIN);
0427       inspect_ordered_towers(hco_towers_from_bestcluster,true_energy,SvtxTrack::HCALOUT);
0428
0429       // ----------------------------------------------------------------------
0430       // ----------------------------------------------------------------------
0431       // ----------------------------------------------------------------------
0432
0433
0434
0435       // --- energy variables from the best candidate clusters
0436       float emc_energy_best = -9999;
0437       float hci_energy_best = -9999;
0438       float hco_energy_best = -9999;
0439
0440       // --- energy variables directly from the track object
0441       float emc_energy_track = -9999;
0442       float hci_energy_track = -9999;
0443       float hco_energy_track = -9999;
0444
0445       if ( verbosity > 2 ) cout << "Now attempting to get the energies..." << endl;
0446
0447       // --- get the energy values directly from the best candidate IF the cluster container exists
0448       if ( emc_bestcluster ) emc_energy_best = emc_bestcluster->get_energy();
0449       if ( hci_bestcluster ) hci_energy_best = hci_bestcluster->get_energy();
0450       if ( hco_bestcluster ) hco_energy_best = hco_bestcluster->get_energy();
0451
0452       if ( verbosity > 5 )
0453     {
0454       cout << "emc_energy_best is " << emc_energy_best << endl;
0455       cout << "hci_energy_best is " << hci_energy_best << endl;
0456       cout << "hco_energy_best is " << hco_energy_best << endl;
0457     }
0458
0459       // --- get the energy values directly from the track IF the cluster container exists
0460       if ( emc_clustercontainer ) emc_energy_track = track->get_cal_cluster_e(SvtxTrack::CEMC);
0461       if ( hci_clustercontainer ) hci_energy_track = track->get_cal_cluster_e(SvtxTrack::HCALIN);
0462       if ( hco_clustercontainer ) hco_energy_track = track->get_cal_cluster_e(SvtxTrack::HCALOUT);
0463
0464       if ( verbosity > 5 )
0465     {
0466       cout << "emc_energy_track is " << emc_energy_track << endl;
0467       cout << "hci_energy_track is " << hci_energy_track << endl;
0468       cout << "hco_energy_track is " << hco_energy_track << endl;
0469     }
0470
0471       // -------------------------------------------------------------------------------------
0472       // --- IMPORTANT NOTE: according to Jin, the clusterizing will gather all towers in the
0473       // ---                 calorimeter, so we need to use the 3x3 instead
0474
0475       if ( emc_clustercontainer ) emc_energy_track = track->get_cal_energy_3x3(SvtxTrack::CEMC);
0476       if ( hci_clustercontainer ) hci_energy_track = track->get_cal_energy_3x3(SvtxTrack::HCALIN);
0477       if ( hco_clustercontainer ) hco_energy_track = track->get_cal_energy_3x3(SvtxTrack::HCALOUT);
0478
0479       if ( verbosity > 0 )
0480     {
0481       cout << "emc_energy_track is " << emc_energy_track << endl;
0482       cout << "hci_energy_track is " << hci_energy_track << endl;
0483       cout << "hco_energy_track is " << hco_energy_track << endl;
0484     }
0485
0486       // -------------------------------------------------------------------------------------
0487       // --- IMPORTANT NOTE: according to Jin, dphi and deta will not work correctly in HIJING
0488
0489       float emc_dphi_track = -9999;
0490       float hci_dphi_track = -9999;
0491       float hco_dphi_track = -9999;
0492       if ( emc_clustercontainer ) emc_dphi_track = track->get_cal_dphi(SvtxTrack::CEMC);
0493       if ( hci_clustercontainer ) hci_dphi_track = track->get_cal_dphi(SvtxTrack::HCALIN);
0494       if ( hco_clustercontainer ) hco_dphi_track = track->get_cal_dphi(SvtxTrack::HCALOUT);
0495
0496       float emc_deta_track = -9999;
0497       float hci_deta_track = -9999;
0498       float hco_deta_track = -9999;
0499       if ( emc_clustercontainer ) emc_deta_track = track->get_cal_deta(SvtxTrack::CEMC);
0500       if ( hci_clustercontainer ) hci_deta_track = track->get_cal_deta(SvtxTrack::HCALIN);
0501       if ( hco_clustercontainer ) hco_deta_track = track->get_cal_deta(SvtxTrack::HCALOUT);
0502
0503       _energy_dphi_emc->Fill(true_energy,emc_dphi_track);
0504       _energy_dphi_hci->Fill(true_energy,hci_dphi_track);
0505       _energy_dphi_hco->Fill(true_energy,hco_dphi_track);
0506
0507       _energy_deta_emc->Fill(true_energy,emc_deta_track);
0508       _energy_deta_hci->Fill(true_energy,hci_deta_track);
0509       _energy_deta_hco->Fill(true_energy,hco_deta_track);
0510
0511       float assoc_dphi = 0.1; // adjust as needed, consider class set method
0512       float assoc_deta = 0.1; // adjust as needed, consider class set method
0513       bool good_emc_assoc = fabs(emc_dphi_track) < assoc_dphi && fabs(emc_deta_track) < assoc_deta;
0514       bool good_hci_assoc = fabs(hci_dphi_track) < assoc_dphi && fabs(hci_deta_track) < assoc_deta;
0515       bool good_hco_assoc = fabs(hco_dphi_track) < assoc_dphi && fabs(hco_deta_track) < assoc_deta;
0516
0517       // ------------------------------------------------------------------------------------------
0518
0519       // --- check the variables
0520       if ( verbosity > 3 )
0521     {
0522       cout << "emc_energy_best is " << emc_energy_best << endl;
0523       cout << "hci_energy_best is " << hci_energy_best << endl;
0524       cout << "hco_energy_best is " << hco_energy_best << endl;
0525       cout << "emc_energy_track is " << emc_energy_track << endl;
0526       cout << "hci_energy_track is " << hci_energy_track << endl;
0527       cout << "hco_energy_track is " << hco_energy_track << endl;
0528       cout << "emc_dphi_track is " << emc_dphi_track << endl;
0529       cout << "hci_dphi_track is " << hci_dphi_track << endl;
0530       cout << "hco_dphi_track is " << hco_dphi_track << endl;
0531       cout << "emc_deta_track is " << emc_deta_track << endl;
0532       cout << "hci_deta_track is " << hci_deta_track << endl;
0533       cout << "hco_deta_track is " << hco_deta_track << endl;
0534     } // check on verbosity
0535
0536       float hct_energy_track = 0;
0537       if ( hci_energy_track >= 0 ) hct_energy_track += hci_energy_track;
0538       if ( hco_energy_track >= 0 ) hct_energy_track += hco_energy_track;
0539
0540       float total_energy_dumb = 0;
0541       if ( emc_energy_track >= 0 ) total_energy_dumb += emc_energy_track;
0542       if ( hci_energy_track >= 0 ) total_energy_dumb += hci_energy_track;
0543       if ( hco_energy_track >= 0 ) total_energy_dumb += hco_energy_track;
0544
0545       float total_energy_smart = 0;
0546       if ( good_emc_assoc ) total_energy_smart += emc_energy_track;
0547       if ( good_hci_assoc ) total_energy_smart += hci_energy_track;
0548       if ( good_hco_assoc ) total_energy_smart += hco_energy_track;
0549
0550
0551
0552       float emc_eratio = emc_energy_track/true_energy;
0553       float hci_eratio = hci_energy_track/true_energy;
0554       float hco_eratio = hco_energy_track/true_energy;
0555       float hct_eratio = hct_energy_track/true_energy;
0556       float tot_dumb_eratio = total_energy_dumb/true_energy;
0557       float tot_smart_eratio = total_energy_smart/true_energy;
0558
0559       _energy_ratio_emc->Fill(true_energy,emc_eratio);
0560       _energy_ratio_hci->Fill(true_energy,hci_eratio);
0561       _energy_ratio_hco->Fill(true_energy,hco_eratio);
0562       _energy_ratio_hct->Fill(true_energy,hct_eratio);
0563       _energy_ratio_tot_dumb->Fill(true_energy,tot_dumb_eratio);
0564       _energy_ratio_tot_smart->Fill(true_energy,tot_smart_eratio);
0565
0566       if ( emc_eratio < 0.1 )
0567     {
0568       _energy_ratio_elb_emc->Fill(true_energy,emc_eratio);
0569       _energy_ratio_elb_hci->Fill(true_energy,hci_eratio);
0570       _energy_ratio_elb_hco->Fill(true_energy,hco_eratio);
0571       _energy_ratio_elb_hct->Fill(true_energy,hct_eratio);
0572       _energy_ratio_elb_tot_dumb->Fill(true_energy,tot_dumb_eratio);
0573       _energy_ratio_elb_tot_smart->Fill(true_energy,tot_smart_eratio);
0574     }
0575
0576       if ( emc_eratio > 0.1 )
0577     {
0578       _energy_ratio_eub_emc->Fill(true_energy,emc_eratio);
0579       _energy_ratio_eub_hci->Fill(true_energy,hci_eratio);
0580       _energy_ratio_eub_hco->Fill(true_energy,hco_eratio);
0581       _energy_ratio_eub_hct->Fill(true_energy,hct_eratio);
0582       _energy_ratio_eub_tot_dumb->Fill(true_energy,tot_dumb_eratio);
0583       _energy_ratio_eub_tot_smart->Fill(true_energy,tot_smart_eratio);
0584     }
0585
0586       // ----------------------------------------------------------------------
0587       // ----------------------------------------------------------------------
0588       // ----------------------------------------------------------------------
0589
0590       //cout << "starting the main part of the truth analysis" << endl;
0591
0592     } // end of loop over truth particles
0593
0594
0595
0596   return Fun4AllReturnCodes::EVENT_OK;
0597
0598 }
0599
0600
0601
0602 int STACalorimeterCharacterization::End(PHCompositeNode *topNode)
0603 {
0604   cout << "End called, " << nevents << " events processed with " << nerrors << " errors and " << nwarnings << " warnings." << endl;
0605   return 0;
0606 }
0607
0608
0609
0610 // --- not sure if it's possible but it'd be nice to do some better function polymorphism here
0611 vector<RawCluster*> STACalorimeterCharacterization::get_ordered_clusters(const RawClusterContainer* clusters)
0612 {
0613
0614   if ( clusters == NULL ) return vector<RawCluster*>();
0615
0616   // getClusters has two options, const and non const, so I use const here
0617   auto range = clusters->getClusters();
0618
0619   map<double,RawCluster*> cluster_map;
0620   for ( auto it = range.first; it != range.second; ++it )
0621     {
0622       RawCluster* cluster = it->second;
0623       double energy = cluster->get_energy();
0624       cluster_map.insert(make_pair(energy,cluster));
0625     }
0626
0627   vector<RawCluster*> cluster_list;
0628   for ( auto rit = cluster_map.rbegin(); rit != cluster_map.rend(); ++rit )
0629     {
0630       RawCluster* cluster = rit->second;
0631       cluster_list.push_back(cluster);
0632     }
0633
0634   return cluster_list;
0635
0636 }
0637
0638
0639
0640 // --- not sure if it's possible but it'd be nice to do some better function polymorphism here
0641 vector<RawTower*> STACalorimeterCharacterization::get_ordered_towers(const RawTowerContainer* towers)
0642 {
0643
0644   if ( towers == NULL ) return vector<RawTower*>();
0645
0646   // getTowers has two options, const and non const, so I use const here
0647   auto range = towers->getTowers();
0648
0649   map<double,RawTower*> tower_map;
0650   for ( auto it = range.first; it != range.second; ++it )
0651     {
0652       RawTower* tower = it->second;
0653       double energy = tower->get_energy();
0654       tower_map.insert(make_pair(energy,tower));
0655     }
0656
0657   vector<RawTower*> tower_list;
0658   for ( auto rit = tower_map.rbegin(); rit != tower_map.rend(); ++rit )
0659     {
0660       RawTower* tower = rit->second;
0661       tower_list.push_back(tower);
0662     }
0663
0664   return tower_list;
0665
0666 }
0667
0668
0669
0670 // --- not sure if it's possible but it'd be nice to do some better function polymorphism here
0671 vector<RawTower*> STACalorimeterCharacterization::get_ordered_towers(RawCluster* cluster, RawTowerContainer* towers)
0672 {
0673
0674   if ( cluster == NULL || towers == NULL ) return vector<RawTower*>();
0675
0676   // note that RawCluster* cannot be const, get_towers is not declared as const function in class header
0677   auto range = cluster->get_towers();
0678
0679   multimap<double,RawTower*> tower_map;
0680   for ( auto it = range.first; it != range.second; ++it )
0681     {
0682       // note that RawTowerContainer cannot be const, getTower is not declared as a const function in class header
0683       RawTower* tower = towers->getTower(it->first);
0684       double energy = tower->get_energy();
0685       tower_map.insert(make_pair(energy,tower));
0686     }
0687
0688   vector<RawTower*> tower_list;
0689   for ( auto rit = tower_map.rbegin(); rit != tower_map.rend(); ++rit )
0690     {
0691       RawTower* tower = rit->second;
0692       tower_list.push_back(tower);
0693     }
0694
0695   return tower_list;
0696
0697 }
0698
0699
0700
0701 void STACalorimeterCharacterization::inspect_ordered_towers(const vector<RawTower*>& towers)
0702 {
0703   inspect_ordered_towers(towers,0,0);
0704 }
0705
0706
0707
0708 void STACalorimeterCharacterization::inspect_ordered_towers(const vector<RawTower*>& towers, int calo_layer)
0709 {
0710   inspect_ordered_towers(towers,0,calo_layer);
0711 }
0712
0713
0714
0715 void STACalorimeterCharacterization::inspect_ordered_towers(const vector<RawTower*>& towers, double true_energy, int calo_layer)
0716 {
0717
0718   if ( towers.size() == 0 ) return;
0719
0720   int nphi = 9999;
0721
0722   if ( calo_layer == SvtxTrack::CEMC ) nphi = 256;
0723   if ( calo_layer == SvtxTrack::HCALIN ) nphi = 64;
0724   if ( calo_layer == SvtxTrack::HCALOUT ) nphi = 64;
0725
0726   RawTower* ctower = towers[0]; // central tower
0727
0728   for ( unsigned int i = 0; i < towers.size(); ++i )
0729     {
0730       RawTower* itower = towers[i];
0731       double ienergy = itower->get_energy();
0732
0733       // if ( verbosity > 5 ) cout << "energy is " << ienergy << " tower address is " << itower << endl;
0734
0735       // get the coordinates
0736       int etabin_center = ctower->get_bineta();
0737       int phibin_center = ctower->get_binphi();
0738       int etabin = itower->get_bineta();
0739       int phibin = itower->get_binphi();
0740       // recenter around central tower
0741       etabin -= etabin_center;
0742       phibin -= phibin_center;
0743       // boundary and periodicity corrections...
0744       if ( phibin > nphi/2 ) phibin -= nphi;
0745       if ( phibin < -nphi/2 ) phibin += nphi;
0746       // if ( verbosity > 6 ) cout << "eta phi coordinates relative to central tower " << etabin << " " << phibin << endl;
0747
0748       // fill some 2d coordinate space histograms
0749       if ( calo_layer == SvtxTrack::CEMC)
0750     {
0751       _towers_3x3_emc->Fill(etabin,phibin,ienergy/true_energy);
0752       _towers_5x5_emc->Fill(etabin,phibin,ienergy/true_energy);
0753       _towers_7x7_emc->Fill(etabin,phibin,ienergy/true_energy);
0754       _towers_9x9_emc->Fill(etabin,phibin,ienergy/true_energy);
0755     }
0756       if ( calo_layer == SvtxTrack::HCALIN)
0757     {
0758       _towers_3x3_hci->Fill(etabin,phibin,ienergy/true_energy);
0759       _towers_5x5_hci->Fill(etabin,phibin,ienergy/true_energy);
0760       _towers_7x7_hci->Fill(etabin,phibin,ienergy/true_energy);
0761       _towers_9x9_hci->Fill(etabin,phibin,ienergy/true_energy);
0762     }
0763       if ( calo_layer == SvtxTrack::HCALOUT)
0764     {
0765       _towers_3x3_hco->Fill(etabin,phibin,ienergy/true_energy);
0766       _towers_5x5_hco->Fill(etabin,phibin,ienergy/true_energy);
0767       _towers_7x7_hco->Fill(etabin,phibin,ienergy/true_energy);
0768       _towers_9x9_hco->Fill(etabin,phibin,ienergy/true_energy);
0769     }
0770
0771     } // loop over towers
0772
0773   // characterize the tower energy with the truth energy
0774   double energy_sumtower[10] = {0};
0775   double energy_singletower = 0;
0776   unsigned int nloop = 10;
0777   if ( towers.size() < nloop ) nloop = towers.size();
0778   for ( unsigned int i = 0; i < nloop; ++i )
0779     {
0780       // get the single and summed tower energies
0781       energy_singletower = towers[i]->get_energy();
0782       energy_sumtower[i] = energy_singletower;
0783       for ( unsigned int j = i; j > 0; --j ) energy_sumtower[j] += energy_sumtower[j-1];
0784
0785       // if ( verbosity > 5 )
0786       //   {
0787       //     cout << "inspecting tower energies" << endl;
0788       //     cout << "single tower energy is " << energy_singletower << endl;
0789       //     cout << "sum tower energy is " << energy_sumtower[i] << endl;
0790       //   }
0791
0792       // fill histograms
0793       if ( calo_layer == SvtxTrack::CEMC )
0794     {
0795       _tower_energy_emc[i]->Fill(true_energy,energy_singletower/true_energy);
0796       _towersum_energy_emc[i]->Fill(true_energy,energy_sumtower[i]/true_energy);
0797     }
0798       if ( calo_layer == SvtxTrack::HCALIN )
0799     {
0800       _tower_energy_hci[i]->Fill(true_energy,energy_singletower/true_energy);
0801       _towersum_energy_hci[i]->Fill(true_energy,energy_sumtower[i]/true_energy);
0802     }
0803       if ( calo_layer == SvtxTrack::HCALOUT )
0804     {
0805       _tower_energy_hco[i]->Fill(true_energy,energy_singletower/true_energy);
0806       _towersum_energy_hco[i]->Fill(true_energy,energy_sumtower[i]/true_energy);
0807     }
0808     } // fixed size loop
0809
0810 } // inspect_ordered_towers