Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2026-04-06 08:09:42 

0001 //____________________________________________________________________________..
0002 //
0003 // This is a template for a Fun4All SubsysReco module with all methods from the
0004 // $OFFLINE_MAIN/include/fun4all/SubsysReco.h baseclass
0005 // You do not have to implement all of them, you can just remove unused methods
0006 // here and in HerwigProductionQAModule.h.
0007 //
0008 // HerwigProductionQAModule(const std::string &name = "HerwigProductionQAModule")
0009 // everything is keyed to HerwigProductionQAModule, duplicate names do work but it makes
0010 // e.g. finding culprits in logs difficult or getting a pointer to the module
0011 // from the command line
0012 //
0013 // HerwigProductionQAModule::~HerwigProductionQAModule()
0014 // this is called when the Fun4AllServer is deleted at the end of running. Be
0015 // mindful what you delete - you do loose ownership of object you put on the node tree
0016 //
0017 // int HerwigProductionQAModule::Init(PHCompositeNode *topNode)
0018 // This method is called when the module is registered with the Fun4AllServer. You
0019 // can create historgrams here or put objects on the node tree but be aware that
0020 // modules which haven't been registered yet did not put antyhing on the node tree
0021 //
0022 // int HerwigProductionQAModule::InitRun(PHCompositeNode *topNode)
0023 // This method is called when the first event is read (or generated). At
0024 // this point the run number is known (which is mainly interesting for raw data
0025 // processing). Also all objects are on the node tree in case your module's action
0026 // depends on what else is around. Last chance to put nodes under the DST Node
0027 // We mix events during readback if branches are added after the first event
0028 //
0029 // int HerwigProductionQAModule::process_event(PHCompositeNode *topNode)
0030 // called for every event. Return codes trigger actions, you find them in
0031 // $OFFLINE_MAIN/include/fun4all/Fun4AllReturnCodes.h
0032 //   everything is good:
0033 //     return Fun4AllReturnCodes::EVENT_OK
0034 //   abort event reconstruction, clear everything and process next event:
0035 //     return Fun4AllReturnCodes::ABORT_EVENT; 
0036 //   proceed but do not save this event in output (needs output manager setting):
0037 //     return Fun4AllReturnCodes::DISCARD_EVENT; 
0038 //   abort processing:
0039 //     return Fun4AllReturnCodes::ABORT_RUN
0040 // all other integers will lead to an error and abort of processing
0041 //
0042 // int HerwigProductionQAModule::ResetEvent(PHCompositeNode *topNode)
0043 // If you have internal data structures (arrays, stl containers) which needs clearing
0044 // after each event, this is the place to do that. The nodes under the DST node are cleared
0045 // by the framework
0046 //
0047 // int HerwigProductionQAModule::EndRun(const int runnumber)
0048 // This method is called at the end of a run when an event from a new run is
0049 // encountered. Useful when analyzing multiple runs (raw data). Also called at
0050 // the end of processing (before the End() method)
0051 //
0052 // int HerwigProductionQAModule::End(PHCompositeNode *topNode)
0053 // This is called at the end of processing. It needs to be called by the macro
0054 // by Fun4AllServer::End(), so do not forget this in your macro
0055 //
0056 // int HerwigProductionQAModule::Reset(PHCompositeNode *topNode)
0057 // not really used - it is called before the dtor is called
0058 //
0059 // void HerwigProductionQAModule::Print(const std::string &what) const
0060 // Called from the command line - useful to print information when you need it
0061 //
0062 //____________________________________________________________________________..
0063 
0064 #include "HerwigProductionQAModule.h"
0065 
0066 #include <fun4all/Fun4AllReturnCodes.h>
0067 
0068 #include <phool/PHCompositeNode.h>
0069 
0070 //____________________________________________________________________________..
0071 HerwigProductionQAModule::HerwigProductionQAModule(const std::string data_type_label, const std::string outfile, float trigger, int verb, const std::string &name):
0072  SubsysReco(name)
0073 , output_file(outfile)
0074 , verbosity(verb)
0075 , trigger_val(trigger)  
0076 {
0077     //A very basic module to check the kinematics of a Herwig production
0078     if( data_type_label.find("Photon") != std::string::npos || data_type_label.find("photon") != std::string::npos || data_type_label.find("MB") != std::string::npos ) 
0079         photon  = true;
0080     if( data_type_label.find("Jet") != std::string::npos || data_type_label.find("jet") != std::string::npos || data_type_label.find("MB") != std::string::npos ) 
0081         jet     = true;
0082     if( data_type_label.find("Herwig") != std::string::npos || data_type_label.find("herwig") != std::string::npos ) 
0083         herwig  = true;
0084     else if ( data_type_label.find("Pythia") != std::string::npos || data_type_label.find("pythia") != std::string::npos) 
0085         pythia  = true;
0086     else 
0087         no_gen  = true; //catch for an issue of having no generator, on first event immediately fail 
0088 }
0089 
0090 //____________________________________________________________________________..
0091 HerwigProductionQAModule::~HerwigProductionQAModule()
0092 {
0093 }
0094 
0095 //____________________________________________________________________________..
0096 int HerwigProductionQAModule::Init(PHCompositeNode *topNode)
0097 {
0098     if(verbosity > 1) std::cout<<"Trigger val: " <<trigger_val <<std::endl;
0099     if(jet)
0100     {
0101         for(int i=2; i<7; i++)
0102         {
0103             //all jets
0104             TH1F* h_all_jet_pt=new TH1F(Form("h_jet_r0%d_pt", i ), 
0105                     Form(" R=0.%d jets; p_{T}^{jet}[GeV]; N_{jets}", i),
0106                     100, -0.5, 99.5);
0107             TH1F* h_all_jet_eta=new TH1F(Form("h_jet_r0%d_eta", i ), 
0108                     Form(" R=0.%d jets; #eta_{jet}; N_{jets}", i),
0109                     100, -1.2, 1.2);
0110             TH1F* h_all_jet_phi=new TH1F(Form("h_jet_r0%d_phi", i ), 
0111                     Form(" R=0.%d jets; #varphi_{jet}; N_{jets}", i),
0112                     100, 0, 2*M_PI);
0113             TH1F* h_all_jet_e=new TH1F(Form("h_jet_r0%d_e", i ), 
0114                     Form(" R=0.%d jets; E_{jet}[GeV]; N_{jets}", i),
0115                     100, -0.5, 99.5);
0116             TH1I* h_all_jet_n_comp = new TH1I(Form("h_jet_r0%d_comp", i),
0117                     Form(" R=0.%d jets; N_{comp}; N_{jets}", i),
0118                     100, 0, 100);
0119             //add to the vector
0120             h_all_jets_pt.push_back(h_all_jet_pt);
0121             h_all_jets_eta.push_back(h_all_jet_eta);
0122             h_all_jets_phi.push_back(h_all_jet_phi);
0123             h_all_jets_e.push_back(h_all_jet_e);
0124             h_all_jets_n_comp.push_back(h_all_jet_n_comp);
0125 
0126             //leading jets
0127             TH1F* h_lead_jet_pt=new TH1F(Form("h_lead_jet_r0%d_pt", i ), 
0128                     Form(" R=0.%d jets; p_{T}^{lead jet}[GeV]; N_{evts}", i),
0129                     100, -0.5, 99.5);
0130             TH1F* h_lead_jet_eta=new TH1F(Form("h_lead_jet_r0%d_eta", i ), 
0131                     Form(" R=0.%d jets; #eta_{jet}^{lead}; N_{evts}", i),
0132                     100, -1.2, 1.2);
0133             TH1F* h_lead_jet_phi=new TH1F(Form("h_lead_jet_r0%d_phi", i ), 
0134                     Form(" R=0.%d jets; #varphi_{jet}^{lead}; N_{evts}", i),
0135                     100, 0, 2*M_PI);
0136             TH1F* h_lead_jet_e=new TH1F(Form("h_lead_jet_r0%d_e", i ), 
0137                     Form(" R=0.%d jets; E_{jet}^{lead}[GeV]; N_{evts}", i),
0138                     100, -0.5, 99.5);
0139             TH1I* h_lead_jet_n_comp = new TH1I(Form("h_lead_jet_r0%d_comp", i),
0140                     Form(" R=0.%d jets; N_{comp}^{lead jet}; N_{evts}", i),
0141                     100, 0, 100);
0142             //add to the vector
0143             h_lead_jets_pt.push_back(h_lead_jet_pt);
0144             h_lead_jets_eta.push_back(h_lead_jet_eta);
0145             h_lead_jets_phi.push_back(h_lead_jet_phi);
0146             h_lead_jets_e.push_back(h_lead_jet_e);
0147             h_lead_jets_n_comp.push_back(h_lead_jet_n_comp);
0148 
0149             TH1I* h_n_jet = new TH1I(Form("h_jet_r0%d_n", i), 
0150                     Form(" R=0.%d jets; N_{jet}; N_{event}", i),
0151                     50, 0, 50);
0152             h_n_jets.push_back(h_n_jet);
0153         }
0154     }
0155     if(photon)
0156     {
0157         //all photons
0158         h_all_photons_pt=new TH1F("h_photons_pt", 
0159                 "all  photons; p_{T}^{photon}[GeV]; N_{photons}",
0160                 100, -0.5, 99.5);
0161         h_all_photons_eta=new TH1F("h_photons_eta", 
0162                 "all  photons; #eta_{photon}; N_{photons}",
0163                 100, -1.2, 1.2);
0164         h_all_photons_phi=new TH1F("h_photons_phi", 
0165                 "all  photons; #varphi_{photon}; N_{photons}",
0166                 100, 0, 2*M_PI);
0167         h_all_photons_e=new TH1F("h_photons_e", 
0168                 "all  photons; E_{photon}[GeV]; N_{photons}",
0169                 100, -0.5, 99.5);
0170         h_n_photons     =new TH1I("h_n_photons",
0171                 "all  photons; N_{photons}; N_{evts}",
0172                 100, 0, 100);
0173         
0174         //leading photons
0175         h_lead_photons_pt=new TH1F("h_lead_photons_pt", 
0176                 "lead  photons; p_{T}^{lead photon}[GeV]; N_{evts}",
0177                 100, -0.5, 99.5);
0178         h_lead_photons_eta=new TH1F("h_lead_photons_eta", 
0179                 "lead  photons; #eta_{photon}^{lead}; N_{evts}",
0180                 100, -1.2, 1.2);
0181         h_lead_photons_phi=new TH1F("h_lead_photons_phi", 
0182                 "lead  photons; #varphi_{photon}^{lead}; N_{evts}",
0183                 100, 0, 2*M_PI);
0184         h_lead_photons_e=new TH1F("h_lead_photons_e", 
0185                 "lead  photons; E_{photon}^{lead}[GeV]; N_{evts}",
0186                 100, -0.5, 99.5);
0187         
0188         //direct photons
0189         h_direct_photons_pt=new TH1F("h_direct_photons_pt", 
0190                 "direct  photons; p_{T}^{photon}[GeV]; N_{photon}",
0191                 100, -0.5, 99.5);
0192         h_direct_photons_eta=new TH1F("h_direct_photons_eta", 
0193                 "direct  photons; #eta_{photon}^{direct}; N_{photon}",
0194                 100, -1.2, 1.2);
0195         h_direct_photons_phi=new TH1F("h_direct_photons_phi", 
0196                 "direct  photons; #varphi_{photon}^{direct}; N_{photon}",
0197                 100, 0, 2*M_PI);
0198         h_direct_photons_e=new TH1F("h_direct_photons_E", 
0199                 "direct  photons; E_{photon}[GeV]; N_{photon}",
0200                 100, -0.5, 99.5);
0201         h_n_direct  =new TH1I("h_n_direct",
0202                 "direct  photons; N_{photons}; N_{photon}",
0203                 100, 0, 100);
0204         
0205         //fragmentation photons
0206         h_frag_photons_pt=new TH1F("h_frag_photons_pt", 
0207                 "fragmentation photons; p_{T}^{photon}[GeV]; N_{photon}",
0208                 100, -0.5, 99.5);
0209         h_frag_photons_eta=new TH1F("h_frag_photons_eta", 
0210                 "fragmentation photons; #eta_{photon}^{frag}; N_{photon}",
0211                 100, -1.2, 1.2);
0212         h_frag_photons_phi=new TH1F("h_frag_photons_phi", 
0213                 "fragmentation photons; #varphi_{photon}^{frag}; N_{photon}",
0214                 100, 0, 2*M_PI);
0215         h_frag_photons_e=new TH1F("h_frag_photons_E", 
0216                 "fragmentation photons; E_{photon}[GeV]; N_{photon}",
0217                 100, -0.5, 99.5);
0218         h_n_frag    =new TH1I("h_n_frag",
0219                 "fragmentation  photons; N_{photons}; N_{photon}",
0220                 100, 0, 100);
0221 
0222     }
0223     if(photon && jet)
0224     {
0225         for(int i=2; i<7; i++)
0226         {
0227             //leading jet + photon
0228             TH2F* h_photon_jet_pt_a=new TH2F(Form("h_photon_jet_r0%d_pt", i ), 
0229                     Form(" R=0.%d jets; p_{T}^{photon}[GeV]; p_{T}^{jet}[GeV]; N_{evts}", i),
0230                     100, -0.5, 99.5, 100, -0.5, 99.5);
0231             TH2F* h_photon_jet_eta_a=new TH2F(Form("h_photon_jet_r0%d_eta", i ), 
0232                     Form(" R=0.%d jets; #eta_{photon}; #eta_{jet}; N_{evts}", i),
0233                     100, -1.2, 1.2, 100, -1.2, 1.2);
0234             TH2F* h_photon_jet_phi_a=new TH2F(Form("h_photon_jet_r0%d_phi", i ), 
0235                     Form(" R=0.%d jets; #varphi_{photon}; #varphi_{jet}; N_{evts}", i),
0236                     100, 0, 2*M_PI, 100, 0, 2*M_PI);
0237             TH2F* h_photon_jet_e_a=new TH2F(Form("h_photon_jet_r0%d_e", i ), 
0238                     Form(" R=0.%d jets; E_{photon}[GeV]; E_{jet}; N_{evts}", i),
0239                     100, -0.5, 99.5, 100, -0.5, 99.5);
0240             TH1F* h_photon_jet_dphi_a = new TH1F(Form("h_photon_jet_r0%d_dphi", i),
0241                     Form(" R=0.%d jets; |#Delta #varphi|_{photon jet}; N_{evts}", i),
0242                     100, -0.1, M_PI+0.1);
0243             //add to the vector
0244             h_photon_jet_pt.push_back(h_photon_jet_pt_a);
0245             h_photon_jet_eta.push_back(h_photon_jet_eta_a);
0246             h_photon_jet_phi.push_back(h_photon_jet_phi_a);
0247             h_photon_jet_e.push_back(h_photon_jet_e_a);
0248             h_photon_jet_dphi.push_back(h_photon_jet_dphi_a);
0249         }
0250     }
0251                 
0252     //event categorization 
0253     //1D distributions
0254     h_particle_eta  = new TH1F("h_particle_eta" , "Final State Particles; #eta; N_{particles}", 100, -1.2, 1.2);
0255     h_particle_phi  = new TH1F("h_particle_phi" , "Final State Particles; #varphi; N_{particles}", 100, 0, 2*M_PI);
0256     h_particle_e    = new TH1F("h_particle_e"   , "Final State Particles; E [GeV]; N_{particles}", 200, -0.5, 99.5);
0257     h_particle_et   = new TH1F("h_particle_et"  , "Final State Particles; E_{T} [GeV]; N_{particles}", 200, -0.5, 99.5);
0258     h_particle_pt   = new TH1F("h_particle_pt"  , "Final State Particles; p_{T} [GeV]; N_{particles}", 200, -0.5, 99.5);
0259     h_total_E   = new TH1F("h_total_E"      , "Final State Particles in |#eta| < 1.1; #sum E [GeV]; N_{evts}", 300, -0.5, 299.5);
0260 
0261     //2D correlations
0262     h_particle_et_eta = new TH2F("h_particle_et_eta", 
0263             "Final State Particles; #eta; E_{T} [GeV]; N_{particles}", 
0264             200, -1.2, 1.2, 100, -0.5, 99.5);
0265     h_particle_et_phi = new TH2F("h_particle_et_phi", 
0266             "Final State Particles; #varphi; E_{T} [GeV]; N_{particles}", 
0267             200, 0, 2*M_PI, 100, -0.5, 99.5);
0268     h_particle_pt_eta = new TH2F("h_particle_pt_eta", 
0269             "Final State Particles; #eta; p_{T} [GeV]; N_{particles}", 
0270             200, -1.2, 1.2, 100, -0.5, 99.5);
0271     h_particle_pt_phi = new TH2F("h_particle_pt_phi", 
0272             "Final State Particles; #varphi; p_{T} [GeV]; N_{particles}", 
0273             200, 0, 2*M_PI, 100, -0.5, 99.5);
0274     h_particle_e_eta = new TH2F("h_particle_e_eta", 
0275             "Final State Particles; #eta; E [GeV]; N_{particles}", 
0276             200, -1.2, 1.2, 100, -0.5, 99.5);
0277     h_particle_e_phi = new TH2F("h_particle_e_phi", 
0278             "Final State Particles; #varphi; E [GeV]; N_{particles}", 
0279             200, 0, 2*M_PI, 100, -0.5, 99.5);
0280     h_particle_phi_eta = new TH2F("h_particle_phi_eta", 
0281             "Final State Particles; #eta; #varphi; N_{particles}", 
0282             100, -1.2, 1.2, 100, 0, 2*M_PI);
0283     
0284     //electrons
0285     h_electron_phi_eta = new TH2F("h_electron_phi_eta", 
0286             "Final State Particles; #eta; #varphi; N_{electrons}", 
0287             100, -1.2, 1.2, 100, 0, 2*M_PI);
0288     h_electron_pt   = new TH1F("h_electron_pt"  , "Final State Particles; p_{T} [GeV]; N_{electrons}", 200, -0.5, 99.5);
0289     
0290     //protons
0291     h_proton_phi_eta = new TH2F("h_proton_phi_eta", 
0292             "Final State Particles; #eta; #varphi; N_{protons}", 
0293             100, -1.2, 1.2, 100, 0, 2*M_PI);
0294     h_proton_pt = new TH1F("h_proton_pt"  , "Final State Particles; p_{T} [GeV]; N_{protons}", 200, -0.5, 99.5);
0295     
0296     //neutrons
0297     h_neutron_phi_eta = new TH2F("h_neutron_phi_eta", 
0298             "Final State Particles; #eta; #varphi; N_{neutrons}", 
0299             100, -1.2, 1.2, 100, 0, 2*M_PI);
0300     h_neutron_pt    = new TH1F("h_neutron_pt"  , "Final State Particles; p_{T} [GeV]; N_{neutrons}", 200, -0.5, 99.5);
0301     
0302     //pions
0303     h_pion_phi_eta  = new TH2F("h_pion_phi_eta", 
0304             "Final State Particles; #eta; #varphi; N_{pions}", 
0305             100, -1.2, 1.2, 100, 0, 2*M_PI);
0306     h_pion_pt   = new TH1F("h_pion_pt"  , "Final State Particles; p_{T} [GeV]; N_{pions}", 200, -0.5, 99.5);
0307         
0308     //photons
0309     h_photon_phi_eta = new TH2F("h_photon_phi_eta", 
0310             "Final State Particles; #eta; #varphi; N_{photons}", 
0311             100, -1.2, 1.2, 100, 0, 2*M_PI);
0312     h_photon_pt = new TH1F("h_photon_pt"  , "Final State Particles; p_{T} [GeV]; N_{photons}", 200, -0.5, 99.5);
0313     //counting events 
0314     h_particle_n    = new TH1I("h_particle_n", "Final state particles; N_{particle}; N_{evts}", 200, 0, 200);
0315     h_electron_n    = new TH1I("h_electron_n", "Final state electrons; N_{electron}; N_{evts}", 200, 0, 200);
0316     h_proton_n  = new TH1I("h_proton_n", "Final state protons; N_{proton}; N_{evts}", 200, 0, 200);
0317     h_neutron_n = new TH1I("h_neutron_n", "Final state neutrons; N_{neutron}; N_{evts}", 200, 0, 200);
0318     h_pion_n    = new TH1I("h_pion_n", "Final state #pi; N_{#pi}; N_{evts}", 200, 0, 200);
0319     h_photon_n  = new TH1I("h_photon_n", "Final state photons; N_{#gamma}; N_{evts}", 200, 0, 200);
0320     
0321     return Fun4AllReturnCodes::EVENT_OK;
0322 }
0323 
0324 //____________________________________________________________________________..
0325 int HerwigProductionQAModule::InitRun(PHCompositeNode *topNode)
0326 {
0327     return Fun4AllReturnCodes::EVENT_OK;
0328 }
0329 
0330 //____________________________________________________________________________..
0331 int HerwigProductionQAModule::process_event(PHCompositeNode *topNode)
0332 {
0333     if(no_gen) return Fun4AllReturnCodes::ABORTRUN;
0334     n_evt++;
0335     if(verbosity >= 1 ) std::cout<<"Working on event " <<n_evt <<std::endl;
0336     if(herwig) process_herwig_event(topNode); //processes just a hepmc node
0337     if(pythia) process_pythia_event(topNode); //processes a fully reconstructed pythia sample
0338     
0339     return Fun4AllReturnCodes::EVENT_OK;
0340 }
0341 int HerwigProductionQAModule::process_herwig_event(PHCompositeNode* topNode){
0342     //process data with a HepMC input 
0343     std::vector<std::vector<Jet*>*> identified_jets {}; //to hold r=0.2-r=0.6 jets 
0344     std::vector<HepMC::GenParticle*> photons {};
0345     std::vector<HepMC::GenParticle*> event_particles {};
0346 //  std::array<float,3> vertex {};
0347     PHHepMCGenEventMap* phg=findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0348     if(!phg){ 
0349         return 1; //catch empty event pmap
0350     }
0351     else{
0352         for(PHHepMCGenEventMap::ConstIter eventIter=phg->begin(); eventIter != phg->end(); ++eventIter){
0353             PHHepMCGenEvent* hepev=eventIter->second;
0354             if(!hepev){
0355                 continue;
0356             }
0357             else{
0358                 HepMC::GenEvent* ev=hepev->getEvent();
0359                 if(!ev) continue;
0360                 else{
0361                     auto vtx = ev->signal_process_vertex();
0362                     if(!vtx) continue;
0363 /*                  else{ //fill in info from the generator vertex
0364                         auto vtx_pos = vtx->position();
0365                         vertex[0] = vtx_pos->x();
0366                         vertex[1] = vtx_pos->y();
0367                         vertex[2] = vtx_pos->z();
0368                         h_vertex_x->Fill(vertex[0]);
0369                         h_vertex_y->Fill(vertex[1]);
0370                         h_vertex_z->Fill(vertex[2]);
0371                         h_vertex_rz->Fill(std::sqrt(std::pow(vertex[0], 2) + std::pow(vertex[1] ,2)), vertex[3]);
0372                         h_vertex_thetaz->Fill(std::atan2(vertex[1], vertex[0]), vertex[2]);
0373                     }*/
0374                     for(HepMC::GenEvent::particle_const_iterator iter=ev->particles_begin(); iter != ev->particles_end(); ++iter) 
0375                     {
0376                         auto particle=(*iter);
0377                         if(!particle) continue;
0378                         if(particle->status()==1 && !(particle->end_vertex()) )
0379                         {
0380                             event_particles.push_back(particle);
0381                             if( particle->pdg_id() == 22)
0382                             {
0383                                 photons.push_back(particle);
0384                             }
0385                         }
0386                     }
0387                 }
0388             }
0389         }
0390     }
0391     findJets(event_particles, &identified_jets);
0392     if(photon)      runAnalysisPhotonJets(identified_jets, photons);
0393     else if (jet)       runAnalysisJets(identified_jets);
0394     if( photon || jet )     runAnalysisEvent(event_particles);  
0395     return Fun4AllReturnCodes::EVENT_OK;
0396 
0397 }
0398 void HerwigProductionQAModule::findJets(std::vector<HepMC::GenParticle*> event_particles, std::vector<std::vector<Jet*>*>* jets)
0399 {
0400     //just invoke fastJet over the relevant event particles and get some sembelance of truth jets at the HepMC level
0401     std::vector<fastjet::PseudoJet> candidates;
0402     fastjet::JetDefinition fjd2 (fastjet::antikt_algorithm, 0.2);
0403     fastjet::JetDefinition fjd3 (fastjet::antikt_algorithm, 0.3);
0404     fastjet::JetDefinition fjd4 (fastjet::antikt_algorithm, 0.4);
0405     fastjet::JetDefinition fjd5 (fastjet::antikt_algorithm, 0.5);
0406     fastjet::JetDefinition fjd6 (fastjet::antikt_algorithm, 0.6);
0407     for(auto p:event_particles)
0408     {
0409         auto pid = p->pdg_id();
0410         if (abs(pid) > 11 &&  abs(pid) < 19 ) continue;
0411         auto mom = p->momentum();
0412         candidates.push_back(fastjet::PseudoJet(mom.px(), mom.py(), mom.pz(), mom.e()));
0413     }   
0414     fastjet::ClusterSequence cs2(candidates, fjd2);
0415     fastjet::ClusterSequence cs3(candidates, fjd3);
0416     fastjet::ClusterSequence cs4(candidates, fjd4);
0417     fastjet::ClusterSequence cs5(candidates, fjd5);
0418     fastjet::ClusterSequence cs6(candidates, fjd6);
0419     auto js2=cs2.inclusive_jets(1.);
0420     auto js3=cs3.inclusive_jets(1.);
0421     auto js4=cs4.inclusive_jets(1.);
0422     auto js5=cs5.inclusive_jets(1.);
0423     auto js6=cs6.inclusive_jets(1.);
0424     std::vector< std::vector<fastjet::PseudoJet>  > jets_size
0425     {
0426         js2,
0427         js3,
0428         js4,
0429         js5,
0430         js6
0431     };
0432     for(auto js:jets_size)
0433     {
0434         JetContainerv1* jc=new JetContainerv1();
0435         std::vector<Jet*>* jt=new std::vector<Jet*>();
0436         for(auto j:js)
0437         {
0438             auto jet = jc->add_jet();
0439             jet->set_px(j.px());
0440             jet->set_py(j.py());
0441             jet->set_pz(j.pz());
0442             jet->set_e( j.e() );
0443             for(auto cmp:j.constituents())
0444             {
0445                 jet->insert_comp(Jet::SRC::HEPMC_IMPORT, cmp.user_index());
0446             }
0447             jt->push_back((Jet*)jet);
0448         }
0449         jets->push_back(jt);
0450     }
0451     return;
0452         
0453 }
0454 
0455 int HerwigProductionQAModule::process_pythia_event(PHCompositeNode* topNode){
0456     //just have to extract the HepMC input part of the DST
0457     //also grab the jets
0458     std::vector<std::vector<Jet*>*>  jets;
0459     std::vector<HepMC::GenParticle*> photons;
0460     std::vector<HepMC::GenParticle*> event_particles;
0461     auto hepmc_gen_event= findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0462     if(!hepmc_gen_event) return Fun4AllReturnCodes::EVENT_OK; 
0463     else if(hepmc_gen_event)
0464     {
0465         for(PHHepMCGenEventMap::ConstIter evtIter=hepmc_gen_event->begin(); evtIter != hepmc_gen_event->end(); ++evtIter)
0466         {
0467             PHHepMCGenEvent* hepev=evtIter->second;
0468             if(!hepev) continue;
0469             else if(hepev){
0470                 HepMC::GenEvent* ev= hepev->getEvent();
0471                 if( !ev ) continue;
0472                 else if( ev )
0473                 {
0474                     for(HepMC::GenEvent::particle_const_iterator iter=ev->particles_begin(); iter != ev->particles_end(); ++iter) 
0475                     {
0476                         auto particle=(*iter);
0477                         if(!particle) continue;
0478                         if(particle->status()==1 && !particle->end_vertex()) //picks up final state particles only 
0479                         {
0480                             event_particles.push_back(particle);
0481                             if( abs(particle->pdg_id()) == 22)
0482                             {
0483                                 photons.push_back(particle);
0484                             }
0485                         }
0486                     }
0487                 }
0488             }
0489         }
0490     }   
0491     auto js2=findNode::getClass<JetContainerv1>(topNode, "AntiKt_Truth_r02");
0492     auto js3=findNode::getClass<JetContainerv1>(topNode, "AntiKt_Truth_r03");
0493     auto js4=findNode::getClass<JetContainerv1>(topNode, "AntiKt_Truth_r04");
0494     auto js5=findNode::getClass<JetContainerv1>(topNode, "AntiKt_Truth_r05");
0495     auto js6=findNode::getClass<JetContainerv1>(topNode, "AntiKt_Truth_r06");
0496     std::vector<JetContainerv1*> jts {js2, js3, js4, js5, js6};
0497     if(photons.size() <= 1) return Fun4AllReturnCodes::EVENT_OK;
0498     for(int i=0; i<(int)jts.size(); i++)
0499     {
0500         auto js=jts.at(i);
0501         if(!js) continue;
0502         else if(js)
0503         {
0504             std::vector<Jet*>* jetsize=new std::vector<Jet*>();
0505             for(auto j:*js)
0506             {
0507                 if(!j) continue;
0508                 else if (j) 
0509                 {
0510                     if(j->get_pt() < 1 ) continue; //apply a 1 gev min pt, just as with herwig jets
0511                     jetsize->push_back(j);
0512                 }
0513             }
0514             jets.push_back(jetsize);
0515         }
0516     }
0517 //  if(photon)      runAnalysisPhotonJets(jets, photons);
0518     /*else*/ if (jet)       runAnalysisJets(jets);
0519     if( photon || jet )     runAnalysisEvent(event_particles);  
0520     std::cout<<__LINE__<<std::endl;
0521     return Fun4AllReturnCodes::EVENT_OK;
0522 }
0523 
0524 std::vector<std::array<float, 4>> HerwigProductionQAModule::runAnalysisJets(std::vector<std::vector<Jet*>*> jets_of_all_sizes)
0525 {
0526     //run analysis of jets 
0527     int i=0;
0528     std::vector<std::array<float, 4>> lead_jet_of_all_sizes {};
0529     for(auto jets:jets_of_all_sizes)
0530     {
0531         float lead_pt   = 0.; 
0532         float lead_eta  = 0.;
0533         float lead_phi  = 0.;
0534         float lead_e    = 0.;
0535         int lead_comp   = 0 ;
0536         for(auto jet:*jets){
0537             if(std::abs(jet->get_eta()) > 1.1 - 0.1*(i+2)) continue; //restrict to acceptance of detector
0538             if(!jet) continue;
0539             float phi = jet->get_phi();
0540             if(phi < 0 ) phi+=2*M_PI;
0541             if(jet->get_pt()  > lead_pt){
0542                     lead_pt     = jet->get_pt();
0543                 lead_eta    = jet->get_eta();
0544                 lead_phi    = phi;
0545                 lead_e      = jet->get_e();
0546                 lead_comp   = (int)((jet->get_comp_vec()).size());
0547             }
0548             h_all_jets_pt.at(i)->   Fill(jet->get_pt());
0549             h_all_jets_eta.at(i)->  Fill(jet->get_eta());
0550             h_all_jets_phi.at(i)->  Fill(phi);
0551             h_all_jets_e.at(i)->    Fill(jet->get_e());
0552             h_all_jets_n_comp.at(i)->Fill((int)((jet->get_comp_vec()).size()));
0553         }
0554         h_lead_jets_pt.at(i)->  Fill(lead_pt);
0555         h_lead_jets_eta.at(i)-> Fill(lead_eta);
0556         h_lead_jets_phi.at(i)-> Fill(lead_phi);
0557         h_lead_jets_e.at(i)->   Fill(lead_e);
0558         h_lead_jets_n_comp.at(i)->Fill(lead_comp);
0559         std::array<float, 4> lead_jet {lead_pt, lead_eta, lead_phi, lead_e};
0560         lead_jet_of_all_sizes.push_back(lead_jet);
0561         h_n_jets.at(i)->Fill((int)jets->size());
0562         i++;
0563     }
0564     return lead_jet_of_all_sizes;
0565 }
0566 
0567 int HerwigProductionQAModule::runAnalysisPhotonJets(std::vector<std::vector<Jet*>*> jets_of_all_sizes, std::vector<HepMC::GenParticle*> photons)
0568 {
0569     //run the analysis of the photons + jets
0570     auto lead_jet_of_all_sizes=runAnalysisJets(jets_of_all_sizes);
0571     float lead_pt   = 0.;
0572     float lead_eta  = 0.;
0573     float lead_phi  = 0.;
0574     float lead_e    = 0.;
0575 
0576     int n_frag  = 0;
0577     int n_direct    = 0;
0578     if(photons.size() <= 1 ) return Fun4AllReturnCodes::EVENT_OK;
0579     for(auto ph:photons)
0580     {
0581         auto p=ph->momentum();
0582         float pt=std::sqrt(std::pow(p.px(), 2)+ std::pow(p.py(), 2));
0583         float phi = p.phi();
0584         if(std::abs(p.pseudoRapidity()) > 1.1) continue; //kinematic region of the detector
0585         if(phi < 0 ) phi += 2*M_PI; 
0586         if(pt > lead_pt) 
0587         {
0588             lead_pt     = pt;
0589             lead_eta    = p.pseudoRapidity();
0590             lead_phi    = phi;
0591             lead_e      = p.e();
0592         }
0593         h_all_photons_pt->  Fill(pt);
0594         h_all_photons_eta-> Fill(p.pseudoRapidity());
0595         h_all_photons_phi-> Fill(phi);
0596         h_all_photons_e->   Fill(p.e());
0597     
0598         HepMC::GenVertex* prod = ph->production_vertex();
0599         HepMC::GenEvent* paren = ph->parent_event();
0600         if(*(paren->signal_process_vertex()) == *prod ){
0601             //this is a direct photon 
0602             n_direct++;
0603             h_direct_photons_pt->   Fill(pt);
0604             h_direct_photons_eta->  Fill(p.pseudoRapidity());
0605             h_direct_photons_phi->  Fill(phi);
0606             h_direct_photons_e->    Fill(p.e());
0607         }
0608         else
0609         {
0610             n_frag++;
0611             h_frag_photons_pt-> Fill(pt);
0612             h_frag_photons_eta->    Fill(p.pseudoRapidity());
0613             h_frag_photons_phi->    Fill(phi);
0614             h_frag_photons_e->  Fill(p.e());
0615         }
0616 
0617     }
0618     h_lead_photons_pt->Fill(lead_pt);
0619     h_lead_photons_eta->Fill(lead_eta);
0620     h_lead_photons_phi->Fill(lead_phi);
0621     h_lead_photons_e->Fill(lead_e);
0622     
0623     h_n_photons->Fill((int)photons.size());
0624     h_n_frag->Fill(n_frag);
0625     h_n_direct->Fill(n_direct);
0626     if(jet){
0627     //now compare correlations between lead jet and lead photon
0628         for(int i=0; i < (int)lead_jet_of_all_sizes.size(); i++)
0629         {
0630                auto lead_jet = lead_jet_of_all_sizes.at(i);
0631                float lead_jet_pt    = lead_jet[0];         
0632                float lead_jet_eta   = lead_jet[1];
0633                float lead_jet_phi   = lead_jet[2];
0634                float lead_jet_e     = lead_jet[3];
0635                h_photon_jet_pt.at(i)->Fill(lead_pt, lead_jet_pt);
0636                h_photon_jet_eta.at(i)->Fill(lead_eta, lead_jet_eta);
0637                h_photon_jet_phi.at(i)->Fill(lead_phi, lead_jet_phi);
0638                h_photon_jet_e.at(i)->Fill(lead_e, lead_jet_e);
0639                float delta_phi = std::abs(lead_phi - lead_jet_phi);
0640                if(delta_phi > M_PI) delta_phi = 2* M_PI - delta_phi; 
0641                h_photon_jet_dphi.at(i)->Fill(delta_phi);
0642         }
0643     }
0644     return Fun4AllReturnCodes::EVENT_OK;
0645 }
0646 int HerwigProductionQAModule::runAnalysisEvent(std::vector<HepMC::GenParticle*> particles)
0647 {
0648     //this is just a QA of bulk properties
0649     float total_E   = 0.;
0650     int n_p     = 0;
0651     int n_e     = 0;
0652     int n_n     = 0;
0653     int n_pi    = 0;
0654     int n_ph    = 0;
0655     h_particle_n->Fill((int)(particles.size()));
0656     for(auto p:particles)
0657     {
0658         float particle_eta  = p->momentum().pseudoRapidity();
0659         float particle_phi  = p->momentum().phi();
0660         float particle_e    = p->momentum().e();
0661         float particle_px   = p->momentum().px();
0662         float particle_py   = p->momentum().pz();
0663         float particle_pt   = std::sqrt(std::pow(particle_px, 2) + std::pow(particle_py, 2));
0664         float particle_et   = particle_e / std::cosh(particle_eta);
0665         int   particle_id   = std::abs(p->pdg_id());
0666         
0667         if(particle_pt < 0.1) continue;
0668         if(std::abs(particle_eta) > 1.1) continue;
0669         if(particle_phi < 0 ) particle_phi += 2*M_PI;
0670         total_E += particle_e;
0671         h_particle_eta->Fill(particle_eta);
0672         h_particle_phi->Fill(particle_phi);
0673         h_particle_e->Fill(particle_e);
0674         h_particle_pt->Fill(particle_pt);
0675         h_particle_et->Fill(particle_et);
0676     
0677         h_particle_et_eta->Fill(particle_eta, particle_et);
0678         h_particle_et_phi->Fill(particle_phi, particle_et);
0679         h_particle_pt_phi->Fill(particle_phi, particle_pt);
0680         h_particle_pt_eta->Fill(particle_eta, particle_pt);
0681         h_particle_phi_eta->Fill(particle_eta, particle_phi);
0682         h_particle_e_phi->Fill(particle_phi, particle_e);
0683         h_particle_e_eta->Fill(particle_eta, particle_e);
0684         if(particle_id == 11 ){
0685                 n_e++;
0686             h_electron_phi_eta->Fill(particle_eta, particle_phi);
0687             h_electron_pt->Fill(particle_pt);
0688         }
0689         else if(particle_id == 111 || particle_id == 211 )
0690         {
0691             n_pi++;
0692             h_pion_phi_eta->Fill(particle_eta, particle_phi);
0693             h_pion_pt->Fill(particle_pt);
0694         }
0695         else if(particle_id == 2212 )
0696         {
0697             n_p++;
0698             h_proton_phi_eta->Fill(particle_eta, particle_phi);
0699             h_proton_pt->Fill(particle_pt);
0700         }
0701         else if(particle_id == 2112)
0702         {
0703             n_n++;
0704             h_neutron_phi_eta->Fill(particle_eta, particle_phi);
0705             h_neutron_pt->Fill(particle_pt);
0706         }
0707         else if(particle_id == 22)
0708         {
0709             n_ph++;
0710             h_photon_phi_eta->Fill(particle_eta, particle_phi);
0711             h_photon_pt->Fill(particle_pt);
0712         }
0713 
0714     }
0715     h_electron_n->Fill(n_e);
0716     h_proton_n->Fill(n_p);
0717     h_neutron_n->Fill(n_n);
0718     h_pion_n->Fill(n_pi);
0719     h_photon_n->Fill(n_p);
0720     h_total_E->Fill(total_E);
0721     return Fun4AllReturnCodes::EVENT_OK;    
0722 }
0723 
0724 //____________________________________________________________________________..
0725 int HerwigProductionQAModule::ResetEvent(PHCompositeNode *topNode)
0726 {
0727     return Fun4AllReturnCodes::EVENT_OK;
0728 }
0729 
0730 //____________________________________________________________________________..
0731 int HerwigProductionQAModule::EndRun(const int runnumber)
0732 {
0733     return Fun4AllReturnCodes::EVENT_OK;
0734 }
0735 
0736 //____________________________________________________________________________..
0737 int HerwigProductionQAModule::End(PHCompositeNode *topNode)
0738 {
0739     return Fun4AllReturnCodes::EVENT_OK;
0740 }
0741 
0742 //____________________________________________________________________________..
0743 int HerwigProductionQAModule::Reset(PHCompositeNode *topNode)
0744 {
0745     return Fun4AllReturnCodes::EVENT_OK;
0746 }
0747 
0748 //____________________________________________________________________________..
0749 void HerwigProductionQAModule::Print(const std::string &what) const
0750 {
0751     TFile* out =new TFile(output_file.c_str(), "RECREATE");
0752     //write out relevant special objects
0753     if(jet)
0754     {
0755         TDirectory* d_j=out->mkdir("Jets");
0756         d_j->cd();
0757         for(int i=0; i<(int) h_all_jets_pt.size(); i++)
0758         {
0759             TDirectory* d_j_i = d_j->mkdir(std::format("R0{}Jets", i+2).c_str(), std::format("R=0.{} Jets", i+2).c_str());
0760             d_j_i->cd();
0761             h_all_jets_pt.at(i)->Write();
0762             h_all_jets_eta.at(i)->Write();
0763             h_all_jets_phi.at(i)->Write();
0764             h_all_jets_e.at(i)->Write();
0765             h_all_jets_n_comp.at(i)->Write();
0766             h_n_jets.at(i)->Write();
0767 
0768             TDirectory* d_j_i_l = d_j_i->mkdir("Lead_Jets", "Lead Jets");   
0769             d_j_i_l->cd();
0770             h_lead_jets_pt.at(i)->Write();
0771             h_lead_jets_eta.at(i)->Write();
0772             h_lead_jets_phi.at(i)->Write();
0773             h_lead_jets_e.at(i)->Write();
0774             h_lead_jets_n_comp.at(i)->Write();
0775     
0776         }
0777     }
0778     if(photon) 
0779     {
0780         TDirectory* d_ph=out->mkdir("Photons");
0781         d_ph->cd();
0782         h_all_photons_pt->Write();  
0783         h_all_photons_eta->Write(); 
0784         h_all_photons_phi->Write(); 
0785         h_all_photons_e->Write();
0786         h_n_photons->Write();
0787 
0788         TDirectory* d_ph_l = d_ph->mkdir("Lead_Photons");   
0789         d_ph_l->cd();
0790         h_lead_photons_pt->Write(); 
0791         h_lead_photons_eta->Write();    
0792         h_lead_photons_phi->Write();    
0793         h_lead_photons_e->Write();
0794         
0795         
0796         TDirectory* d_ph_d=d_ph->mkdir("Direct_Photons");
0797         d_ph_d->cd();   
0798         h_direct_photons_pt->Write();   
0799         h_direct_photons_eta->Write();  
0800         h_direct_photons_phi->Write();  
0801         h_direct_photons_e->Write();
0802         
0803         TDirectory* d_ph_f=d_ph->mkdir("Fragmentation_Photons");
0804         d_ph_f->cd();   
0805         h_frag_photons_pt->Write(); 
0806         h_frag_photons_eta->Write();    
0807         h_frag_photons_phi->Write();    
0808         h_frag_photons_e->Write();
0809 
0810         h_n_frag->Write();
0811         h_n_direct->Write();
0812     }
0813     if(photon && jet)
0814     {
0815         TDirectory* d_phj=out->mkdir("Photon-Jets");
0816         d_phj->cd();
0817         for(int i=0; i<(int)h_photon_jet_pt.size(); i++)
0818         {
0819             TDirectory* d_j_i = d_phj->mkdir(std::format("R0{}Jets", i+2).c_str(), std::format("R=0.{} Jets", i+2).c_str());
0820             d_j_i->cd();
0821             h_photon_jet_pt.at(i)->Write();
0822             h_photon_jet_eta.at(i)->Write();
0823             h_photon_jet_phi.at(i)->Write();
0824             h_photon_jet_e.at(i)->Write();
0825             h_photon_jet_dphi.at(i)->Write();   
0826         }
0827     }
0828     //write out relevant event data 
0829     TDirectory* d_event=out->mkdir("Event");
0830     d_event->cd();
0831     
0832     h_particle_eta->Write();
0833     h_particle_phi->Write();
0834     h_particle_e->Write();
0835     h_particle_pt->Write();
0836     h_particle_et->Write();
0837     h_particle_n->Write();
0838     h_total_E->Write();
0839     
0840     h_particle_et_eta->Write();
0841     h_particle_et_phi->Write();
0842     h_particle_pt_eta->Write();
0843     h_particle_pt_phi->Write();
0844     h_particle_e_phi->Write();
0845     h_particle_e_eta->Write();
0846     h_particle_phi_eta->Write();
0847     
0848     TDirectory* d_part=d_event->mkdir("ParticleCat", "Individual particles");
0849     d_part->cd();
0850     
0851     h_electron_phi_eta->Write();
0852     h_electron_pt->Write();
0853     h_electron_n->Write();
0854         
0855     h_proton_phi_eta->Write();
0856     h_proton_pt->Write();
0857     h_proton_n->Write();
0858 
0859     h_neutron_phi_eta->Write();
0860     h_neutron_pt->Write();
0861     h_neutron_n->Write();
0862 
0863     h_pion_phi_eta->Write();
0864     h_pion_pt->Write();
0865     h_pion_n->Write();
0866     
0867     h_photon_phi_eta->Write();
0868     h_photon_pt->Write();
0869     h_photon_n->Write();
0870 
0871     out->cd();
0872     out->Write();
0873     out->Close();
0874     return;
0875 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!