analysis/TreeMaker/TreeMakerCopyAndMakeJets.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 // --- jet specific stuff
0013 #include <g4jets/Jet.h>
0014 #include <g4jets/JetMap.h>
0015 #include <g4jets/JetV1.h>
0016 #include <g4jets/JetMapV1.h>
0017
0018 using std::cout;
0019 using std::endl;
0020
0021
0022
0023 int TreeMaker::CopyAndMakeJets(PHCompositeNode *topNode)
0024 {
0025
0026   // -----------------------------------------------------------------------------------------------------
0027   // ---
0028   // -----------------------------------------------------------------------------------------------------
0029
0030   // --- calorimeter tower containers
0031   RawTowerContainer* towersEM3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
0032   RawTowerContainer* towersIH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
0033   RawTowerContainer* towersOH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
0034
0035   // --- calorimeter geometry objects
0036   RawTowerGeomContainer* geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0037   RawTowerGeomContainer* geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0038   RawTowerGeomContainer* geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0039
0040   // --- jet objects
0041   JetMap* old_jets = findNode::getClass<JetMap>(topNode,"AntiKt_Tower_r02");
0042   JetMap* new_jets = findNode::getClass<JetMap>(topNode,"AntiKt_Tower_Mod_r02"); // this node is created in createnode
0043   if ( verbosity > 0 )
0044     {
0045       cout << "Regular jet node: " << old_jets << endl;
0046       cout << "Modified jet node: " << new_jets << endl;
0047     }
0048   if ( !old_jets || ! new_jets )
0049     {
0050       cout << "One or more invalid pointers, exiting event" << endl;
0051       return 0;
0052     }
0053
0054   // --- print sizes of old and new objects for diagnostic purposes
0055   if ( verbosity > 0 )
0056     {
0057       cout << "process_event: entering with # original jets = " << old_jets->size() << endl;
0058       cout << "process_event: entering with # new jets = " << new_jets->size() << endl;
0059     }
0060
0061   // --- iterate over jets
0062   int ijet = 0;
0063   for ( JetMap::Iter iter = old_jets->begin(); iter != old_jets->end(); ++iter )
0064     {
0065       // --- get information from this jet (an element of the jet container)
0066       Jet* this_jet = iter->second;
0067       float this_e = this_jet->get_e();
0068       float this_pt = this_jet->get_pt();
0069       float this_phi = this_jet->get_phi();
0070       float this_eta = this_jet->get_eta();
0071       // --- create a new jet object
0072       Jet *new_jet = new JetV1();
0073       float new_total_px = 0;
0074       float new_total_py = 0;
0075       float new_total_pz = 0;
0076       float new_total_e = 0;
0077       if ( verbosity > 3 && this_pt > 5 )
0078         {
0079           cout << "process_event: unsubtracted jet with e / pt / eta / phi = "
0080                << this_e << " / " << this_pt << " / " << this_eta << " / " << this_phi << endl;
0081         }
0082       // --- now loop over individual constituents of this jet
0083       for ( Jet::ConstIter comp = this_jet->begin_comp(); comp != this_jet->end_comp(); ++comp )
0084         {
0085           // --- (*.comp).first is the layer:
0086           // ---  0 = void
0087           // ---  1 = particle
0088           // ---  2 = track
0089           // ---  3 = emcal tower
0090           // ---  4 = emcal cluster
0091           // ---  5 = ihcal tower
0092           // ---  6 = ihcal cluster
0093           // ---  7 = ohcal tower
0094           // ---  8 = ohcal cluster
0095           // ---  9 = femc tower
0096           // --- 10 = femc cluster
0097           // --- 11 = fhcal tower
0098           // --- 12 = fhcal cluster
0099           // --- 13 = emcal retower (combined towers to match ihcal geometry)
0100           // --- 14 = subtracted emcal tower
0101           // --- 15 = subtracted ihcal tower
0102           // --- 16 = subtracted ohcal tower
0103           RawTower *tower = 0;
0104           RawTowerGeom *tower_geom = 0;
0105           double comp_e = 0;
0106           double comp_eta = 0;
0107           double comp_phi = 0;
0108           int comp_ieta = 0;
0109           int comp_iphi = 0;
0110           //double comp_background = 0;
0111           // --- layer 5 is inner hcal towers
0112           if ( (*comp).first == 5 )
0113             {
0114               tower = towersIH3->getTower( (*comp).second );
0115               tower_geom = geomIH->get_tower_geometry(tower->get_key());
0116               comp_ieta = geomIH->get_etabin( tower_geom->get_eta() );
0117               comp_iphi = geomIH->get_phibin( tower_geom->get_phi() );
0118             }
0119           // --- layer 7 is outer hcal towers
0120           else if ( (*comp).first == 7 )
0121             {
0122               tower = towersOH3->getTower( (*comp).second );
0123               tower_geom = geomOH->get_tower_geometry(tower->get_key());
0124               comp_ieta = geomOH->get_etabin( tower_geom->get_eta() );
0125               comp_iphi = geomOH->get_phibin( tower_geom->get_phi() );
0126             }
0127           // --- layer 3 is EMCal towers
0128           // --- layer 13 is re-towered emcal, which uses IHCal geom... not using for now
0129           else if ( (*comp).first == 3 )
0130             {
0131               tower = towersEM3->getTower( (*comp).second );
0132               tower_geom = geomEM->get_tower_geometry(tower->get_key());
0133               comp_ieta = geomEM->get_etabin( tower_geom->get_eta() );
0134               comp_iphi = geomEM->get_phibin( tower_geom->get_phi() );
0135             }
0136
0137           // --- if tower and tower_geom exist, get energy from there
0138           if ( tower && tower_geom )
0139             {
0140               comp_e = tower->get_energy();
0141               comp_eta = tower_geom->get_eta();
0142               comp_phi = tower_geom->get_phi();
0143             }
0144
0145           // --- if very high verbosity, print lots of stuff to screen
0146           if ( verbosity > 4 && this_jet->get_pt() > 5 )
0147             {
0148               cout << "process_event: --> constituent in layer " << (*comp).first
0149                    << ", has unsub E = " << comp_e << ", is at ieta #" << comp_ieta
0150                    << " and iphi # = " << comp_iphi << endl;
0151             }
0152
0153           // --- update constituent energy based on some algorithm to be developed
0154           double comp_new_e = comp_e;
0155
0156           // --- define new kinematics for constituent
0157           double comp_px = comp_new_e / cosh( comp_eta ) * cos( comp_phi );
0158           double comp_py = comp_new_e / cosh( comp_eta ) * sin( comp_phi );
0159           double comp_pz = comp_new_e * tanh( comp_eta );
0160
0161           // --- add up the total for the new jet based on the modified constituents
0162           new_total_px += comp_px;
0163           new_total_py += comp_py;
0164           new_total_pz += comp_pz;
0165           new_total_e += comp_new_e;
0166         } // end of loop over jet constituents
0167
0168       // --- set the properties for the new jet
0169       new_jet->set_px( new_total_px );
0170       new_jet->set_py( new_total_py );
0171       new_jet->set_pz( new_total_pz );
0172       new_jet->set_e( new_total_e );
0173       new_jet->set_id(ijet);
0174
0175       // --- put the jet into a new container class
0176       new_jets->insert( new_jet );
0177
0178       if (verbosity > 1 && this_pt > 5)
0179         {
0180           cout << "process_event: old jet #" << ijet << ", old px / py / pz / e = "
0181                << this_jet->get_px() << " / " << this_jet->get_py() << " / "
0182                << this_jet->get_pz() << " / " << this_jet->get_e() << endl;
0183           cout << "process_event: new jet #" << ijet << ", new px / py / pz / e = "
0184                << new_jet->get_px() << " / " << new_jet->get_py() << " / "
0185                << new_jet->get_pz() << " / " << new_jet->get_e() << endl;
0186         }
0187
0188       ijet++; // increment jet counter
0189
0190     } // end of loop over jet maps (list of jets)
0191
0192   if ( verbosity > 0 )
0193     {
0194       cout << "process_event: exiting with # original jets = " << old_jets->size() << endl;
0195       cout << "process_event: exiting with # new jets = " << new_jets->size() << endl;
0196     }
0197
0198   return 0;
0199
0200 }