Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:13:13 

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 ThirdJetSpectra.h.
0007 //
0008 // ThirdJetSpectra(const std::string &name = "ThirdJetSpectra")
0009 // everything is keyed to ThirdJetSpectra, 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 // ThirdJetSpectra::~ThirdJetSpectra()
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 ThirdJetSpectra::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 ThirdJetSpectra::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 ThirdJetSpectra::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 ThirdJetSpectra::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 ThirdJetSpectra::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 ThirdJetSpectra::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 ThirdJetSpectra::Reset(PHCompositeNode *topNode)
0057 // not really used - it is called before the dtor is called
0058 //
0059 // void ThirdJetSpectra::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 "ThirdJetSpectra.h"
0065 
0066 #include <fun4all/Fun4AllReturnCodes.h>
0067 
0068 #include <phool/PHCompositeNode.h>
0069 
0070 //____________________________________________________________________________..
0071 ThirdJetSpectra::ThirdJetSpectra(std::string numb, const std::string &name):
0072  SubsysReco(name)
0073 {
0074         this->segment_numb=numb;
0075   std::cout << "ThirdJetSpectra::ThirdJetSpectra(const std::string &name) Calling ctor" << std::endl;
0076         xj = new TH1F("xj", "Dijet Momentum Imballance leading to subleading; x_{j}; N_{counts}", 20, 0,1);
0077         xj_strict =new TH1F("xj_strict", "Dijet Momentum imballance leading to subleading (Events passing Dijet Event Cuts); x_{j}; N_{counts}", 20, 0, 1);
0078         xj_onl =new TH1F("xj_onl", "Dijet Momentum imballance leading to subleading without a third jet; x_{j}; N_{counts}", 20, 0, 1);
0079         xj_strict_onl =new TH1F("xj_strict_onl", "Dijet Momentum imballance leading to subleading without a third jet (Events passing Dijet Event Cuts); x_{j}; N_{counts}", 20, 0, 1);
0080         first_jet_pt=new TH1F("first_jet_pt", "p_{T} of leading jet; p_{T} [GeV/c]; N_{counts}", 60, -0.5, 59.5);
0081         second_jet_pt=new TH1F("second_jet_pt", "p_{T} of subleading jet; p_{T} [GeV/c]; N_{counts}", 60, -0.5, 59.5);
0082         third_jet_pt=new TH1F("third_jet_pt", "p_{T} of subsubleading jet; p_{T} [GeV/c]; N_{counts}", 60, -0.5, 59.5);
0083         first_jet_E=new TH1F("first_jet_E", "E of leading jet; E [GeV/c]; N_{counts}", 60, -0.5, 59.5);
0084         second_jet_E=new TH1F("second_jet_E", "E of subleading jet; E [GeV/c]; N_{counts}", 60, -0.5, 59.5);
0085         third_jet_E=new TH1F("third_jet_E", "E of subsubleading jet; E [GeV/c]; N_{counts}", 60, -0.5, 59.5);
0086         first_jet_phi=new TH1F("first_jet_phi", "#varphi of leading jet; #varphi; N_{counts}", 64, -0.5, 2*PI);
0087         second_jet_phi=new TH1F("second_jet_phi", "#varphi of leading jet; #varphi; N_{counts}", 64, 0, 2*PI);
0088         third_jet_phi=new TH1F("third_jet_phi", "#varphi of leading jet; #varphi; N_{counts}", 64, 0, 2*PI);
0089         first_jet_eta=new TH1F("first_jet_eta", "#eta of leading jet; #eta; N_{counts}", 24, -1.11, 1.11);
0090         second_jet_eta=new TH1F("second_jet_eta", "#eta of leading jet; #eta; N_{counts}", 24, -1.11, 1.11);
0091         third_jet_eta=new TH1F("third_jet_eta", "#eta of leading jet; #eta; N_{counts}", 24, -1.11, 1.11);
0092         dphi_12=new TH1F("dphi_12", "#Delta #varphi of leading jet to subleading; #Delta #varphi; N_{counts}", 64, -PI, PI);
0093         dphi_13=new TH1F("dphi_13", "#Delta #varphi of leading jet to subsubleading; #Delta #varphi; N_{counts}", 64, -PI, PI);
0094         dphi_23=new TH1F("dphi_23", "#Delta #varphi of subleading jet to subsubleading;  #Delta #varphi; N_{counts}", 64, -PI, PI);
0095         xj_12 =new TH1F("xj_12", "Dijet Momentum imballance leading to subleading (with third); x_{j}; N_{counts}", 100, 0, 1);
0096         xj_13 =new TH1F("xj_13", "Dijet Momentum imballance leading to subsubleading; x_{j}; N_{counts}", 100, 0, 1);
0097         xj_23 =new TH1F("xj_23", "Dijet Momentum imballance subleading to subsubleading; x_{j}; N_{counts}", 100, 0, 1);
0098         xj_strict_12 =new TH1F("xj_strict_12", "Dijet Momentum imballance leading to subleading with third (Events passing Dijet Event Cuts); x_{j}; N_{counts}", 100, 0, 1);
0099         xj_strict_13 =new TH1F("xj_strict_13", "Dijet Momentum imballance leading to subsubleading (Events passing Dijet Event Cuts); x_{j}; N_{counts}", 100, 0, 1);
0100         xj_strict_23 =new TH1F("xj_strict_23", "Dijet Momentum imballance subleading to subsubleading (Events passing Dijet Event Cuts); x_{j}; N_{counts}", 100, 0, 1);
0101         third_jet_pt_dec = new TH1F("third_jet_pt_dec", "Average Decorrelation of leading and subleading jet as a function of subsubleading jet p_{T}; p_{T} [GeV]; < |#pi - #Delta #varphi| >", 60, -0.5, 59.5);   
0102         third_jet_pt_dec_strict = new TH1F("third_jet_pt_dec_strict", "Average Decorrelation of leading and subleading jet as a function of subsubleading jet p_{T} (Events passing Dijet Event Cuts); p_{T} [GeV]; < |#pi - #Delta #varphi| >", 60, -0.5, 59.5);   
0103         decorr = new TH1F("decorr", "Decorrelation of leading to subleading jets; | #pi -#Delta #varphi |; N_{counts}", 64, 0, PI); 
0104         decorr_strict = new TH1F("decorr_strict", "Decorrelation of leading to subleading jets (Events passing Dijet Event Cuts); | #pi -#Delta #varphi |; N_{counts}", 64, 0, PI);
0105         third_jet_pt_dec_n = new TH2F( "third_jet_pt_dec_n", "Decorrelation of leadidng to subleading jets as fuction of third jet pt; p_{T}^{ssl-jet} [GeV]; |#pi - #Delta #varphi |; N_{counts}", 60, -0.5, 59.5, 64, 0, PI);
0106         third_jet_pt_dec_strict_n = new TH2F("third_jet_pt_dec_strict_n", "Decorrelation of leadidng to subleading jets as fuction of third jet pt (Events passing Dijet Event Cuts); p_{T}^{ssl-jet} [GeV]; |#pi - #Delta #varphi |; N_{counts}", 60, -0.5, 59.5,64,  0, PI);
0107 
0108         dphi_123=new TH3F( "dphi_123", "azimuthal angle of jet objects; #varphi_{leading}; #varphi_{subleading}; #varphi_{subsubleading}; N_{events}", 64, 0, 2*PI, 64, 0, 2*PI, 64, 0, 2*PI);
0109         pt_123= new TH3F("pt123", "Transverse momentum of the jets; p_{T}^{leading} [GeV/c]; p_{T}^{subleading}[GeV/c]; p_{T}^{subsubleading} [GeV/c]; N_{events}", 60, -0.5, 59.5, 60, -0.5, 59.5, 60, -0.5, 59.5); 
0110         DiJEC=new DijetEventCuts(30., 10., 1000., PI/2. );  
0111         DiJEC_strict=new DijetEventCuts(15., 8., 0.7, 3*PI/4. ); //actual kinematics    
0112 }
0113 
0114 //____________________________________________________________________________..
0115 ThirdJetSpectra::~ThirdJetSpectra()
0116 {
0117   std::cout << "ThirdJetSpectra::~ThirdJetSpectra() Calling dtor" << std::endl;
0118 }
0119 
0120 //____________________________________________________________________________..
0121 int ThirdJetSpectra::Init(PHCompositeNode *topNode)
0122 {
0123   std::cout << "ThirdJetSpectra::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0124   return Fun4AllReturnCodes::EVENT_OK;
0125 }
0126 
0127 //____________________________________________________________________________..
0128 int ThirdJetSpectra::InitRun(PHCompositeNode *topNode)
0129 {
0130   std::cout << "ThirdJetSpectra::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0131   return Fun4AllReturnCodes::EVENT_OK;
0132 }
0133 
0134 std::vector<fastjet::PseudoJet> ThirdJetSpectra::findAllJets(HepMC::GenEvent* e1)
0135 {
0136     //do the fast jet clustering, antikt r=-0.4
0137     fastjet::JetDefinition jetdef(fastjet::antikt_algorithm, 0.4);
0138     std::vector<fastjet::PseudoJet> input, output;
0139     for(HepMC::GenEvent::particle_const_iterator iter = e1->particles_begin(); iter != e1->particles_end(); ++iter)
0140     {
0141         if(!(*iter)->end_vertex() && (*iter)->status() == 1){
0142             auto p=(*iter)->momentum();
0143             fastjet::PseudoJet pj( p.px(), p.py(), p.pz(), p.e());
0144             pj.set_user_index((*iter)->barcode());
0145             input.push_back(pj);
0146         }
0147     }
0148     if(input.size() == 0 ) return input;
0149     fastjet::ClusterSequence js (input, jetdef);
0150     auto j = fastjet::sorted_by_pt(js.inclusive_jets(3.));
0151     for(auto j1:j){
0152         if(j1.pt() < 3.) continue;
0153         /*if(j1.pt() < output.back.pt())*/ output.push_back(j1); //just keep in the corect format
0154     }
0155             
0156     return output;
0157 }
0158 void ThirdJetSpectra::getJetTripplet(JetContainerv1* jc, std::vector<Jetv2*>* ds, bool strict)
0159 {
0160     float dphi=PI/2.;
0161     if(strict){
0162          dphi=3*PI/4.;
0163     }
0164     
0165     for(auto j1:*jc){
0166         for(auto j2:*jc){
0167             if(j2->get_pt() >= j1->get_pt() ) continue;
0168             float p1=j1->get_phi();
0169             float p2=j2->get_phi();
0170             if(std::abs(p1-p2) > dphi && std::abs(p1-p2) < PI+dphi ){
0171                 ds->push_back((Jetv2*)j1);
0172                 ds->push_back((Jetv2*)j2);
0173             }
0174             float maxpt=0;
0175             Jetv2* thj;
0176             for(auto j3:*jc){
0177                 if(j3->get_phi() == p1 || j3->get_phi() == p2) continue;
0178                 if(j3->get_pt() > maxpt){
0179                     maxpt=j3->get_pt();
0180                     thj=(Jetv2*)j3;
0181                 }
0182             }
0183             ds->push_back(thj);
0184         }
0185     }
0186     return;
0187 }
0188 void ThirdJetSpectra::getJetPair(JetContainerv1* jc, std::vector<Jetv2*>* ds, bool strict) 
0189 { 
0190     float dphi=PI/2.;
0191     if(strict){
0192         dphi=3*PI/4.;
0193     }
0194     for(auto j1:*jc){
0195         for(auto j2:*jc){
0196             float p1=j1->get_phi();
0197             float p2=j2->get_phi();
0198             if(j2->get_pt() >= j1->get_pt() ) continue;
0199             if(std::abs(p1-p2) > dphi && std::abs(p1-p2) < PI+dphi ){
0200                 ds->push_back((Jetv2*)j1);
0201                 ds->push_back((Jetv2*)j2);
0202             }
0203         }
0204     }
0205     return;
0206 }
0207 //____________________________________________________________________________..
0208 int ThirdJetSpectra::process_event(PHCompositeNode *topNode)
0209 {
0210     PHHepMCGenEventMap* em = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0211     if(!em) return Fun4AllReturnCodes::EVENT_OK;
0212 
0213     for(PHHepMCGenEventMap::ConstIter eventIter = em->begin(); eventIter != em->end(); ++eventIter)
0214     {
0215         bool is_three=false;
0216         PHHepMCGenEvent* pHe=eventIter->second;
0217         if(!pHe ||  pHe->get_embedding_id() != 0 ) continue;
0218         HepMC::GenEvent* truthEvent = pHe->getEvent();
0219         if(!truthEvent) continue;
0220         auto jets=findAllJets(truthEvent);
0221         if(jets.size() < 2) continue;
0222         n_try++;
0223         std::cout<<jets.at(0).pt() <<" is lead jet pt" <<std::endl;
0224         JetContainerv1* jc=new JetContainerv1();
0225         std::vector<float> fake_ratio; //this is just to be able to use my dijet cut object
0226         for(auto i:jets){
0227             fake_ratio.push_back(0.5);  
0228             auto j = jc->add_jet();
0229             j->set_px(i.px());
0230             j->set_py(i.py());
0231             j->set_pz(i.pz());
0232             j->set_e(i.e());
0233             //for(auto cmp: i.constituents()){j->insert_comp(Jet::SRC::PARTICLE, cmp.user_index());}
0234         }
0235         std::array<float, 3> vertex {0,0,0};    
0236         bool is_good_loose=DiJEC->passesTheCut(jc, fake_ratio, fake_ratio, 0.5, vertex, fake_ratio);
0237         bool is_good_strict=DiJEC_strict->passesTheCut(jc, fake_ratio, fake_ratio, 0.5, vertex, fake_ratio);    
0238         if(jets.size() >= 3) pt_123->Fill(jets.at(0).pt(), jets.at(1).pt(), jets.at(2).pt());
0239         if(!is_good_loose) continue;
0240         if (jets.size() >= 3){
0241             n_three++;
0242             is_three=true;
0243         }
0244         n_evts++;
0245         std::vector<Jetv2*> dijet_set, dijet_set_strict;
0246         if(is_three){
0247             getJetTripplet(jc, &dijet_set, false);
0248             if(is_good_strict) getJetTripplet(jc, &dijet_set_strict, true);
0249         }
0250         else{
0251             getJetPair(jc, &dijet_set, false);
0252             if(is_good_strict) getJetPair(jc, &dijet_set_strict, true);
0253         }
0254         float pt1=dijet_set[0]->get_pt();
0255         float pt2=dijet_set[1]->get_pt();
0256         Jetv2* j1=dijet_set[0], *j2=dijet_set[1];
0257         float phi1=j1->get_phi(), phi2=j2->get_phi();
0258         first_jet_pt->Fill(pt1);
0259         first_jet_phi->Fill(j1->get_phi());
0260         first_jet_eta->Fill(j1->get_eta());
0261         first_jet_E->Fill(j1->get_e());
0262         second_jet_pt->Fill(pt2);
0263         second_jet_phi->Fill(j2->get_phi());
0264         second_jet_eta->Fill(j2->get_eta());
0265         second_jet_E->Fill(j2->get_e());
0266         dphi_12->Fill(phi1 - phi2);
0267         float x_j_12 = pt2/pt1;
0268         xj->Fill(x_j_12);
0269         if(!is_three) xj_onl->Fill(x_j_12);
0270         if(is_good_strict){
0271             float xjs_12=dijet_set_strict[1]->get_pt()/dijet_set_strict[0]->get_pt();
0272             xj_strict_12->Fill(xjs_12);
0273             xj_strict->Fill(xjs_12);
0274             if(!is_three) xj_strict_onl->Fill(xjs_12);
0275         }
0276         if(is_three){
0277             Jetv2* j3=dijet_set[3];
0278             float pt3=j3->get_pt();
0279             float phi3=j3->get_phi();
0280             std::cout<<x_j_12<<std::endl;
0281             xj_12->Fill(x_j_12);
0282             third_jet_pt->Fill(pt3);
0283             third_jet_phi->Fill(phi3);
0284             third_jet_eta->Fill(j3->get_eta());
0285             third_jet_E->Fill(j3->get_e());
0286             dphi_13->Fill(phi1-phi3);
0287             dphi_23->Fill(phi2-phi3);
0288             float x13=pt3/pt1;
0289             float x23=pt3/pt2;
0290             xj_13->Fill(x13);
0291             xj_23->Fill(x23);
0292             float decor=std::abs(PI - phi1+phi2);
0293             decorr->Fill(decor);
0294             dphi_123->Fill(std::abs(phi1-phi2), std::abs(phi1-phi3), std::abs(phi2-phi3));
0295             third_jet_pt_dec_n->Fill(pt3, decor);
0296             third_jet_pt_dec->Fill(pt3, decor);
0297             if(is_good_strict){
0298                 Jetv2* j1_s=dijet_set_strict[0];
0299                 Jetv2* j2_s=dijet_set_strict[1];
0300                 Jetv2* j3_s=dijet_set_strict[2];
0301                 float xjs_13=j3_s->get_pt()/j1_s->get_pt();
0302                 float xjs_23=j3_s->get_pt()/j2_s->get_pt();
0303                 float dc_s=std::abs(PI-j1_s->get_phi()+j2_s->get_phi());
0304                 xj_strict_13->Fill(xjs_13);
0305                 xj_strict_23->Fill(xjs_23);
0306                 third_jet_pt_dec_strict->Fill(j3_s->get_pt(), dc_s);    
0307                 third_jet_pt_dec_strict_n->Fill(j3_s->get_pt(), dc_s);  
0308                 decorr_strict->Fill(dc_s);
0309             }
0310         }
0311         
0312     }
0313     return Fun4AllReturnCodes::EVENT_OK;
0314 }
0315 
0316 //____________________________________________________________________________..
0317 int ThirdJetSpectra::ResetEvent(PHCompositeNode *topNode)
0318 {
0319   //std::cout << "ThirdJetSpectra::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0320   return Fun4AllReturnCodes::EVENT_OK;
0321 }
0322 
0323 //____________________________________________________________________________..
0324 int ThirdJetSpectra::EndRun(const int runnumber)
0325 {
0326   std::cout << "ThirdJetSpectra::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0327   return Fun4AllReturnCodes::EVENT_OK;
0328 }
0329 
0330 //____________________________________________________________________________..
0331 int ThirdJetSpectra::End(PHCompositeNode *topNode)
0332 {
0333   std::cout << "ThirdJetSpectra::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0334   return Fun4AllReturnCodes::EVENT_OK;
0335 }
0336 
0337 //____________________________________________________________________________..
0338 int ThirdJetSpectra::Reset(PHCompositeNode *topNode)
0339 {
0340 // std::cout << "ThirdJetSpectra::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0341   return Fun4AllReturnCodes::EVENT_OK;
0342 }
0343 
0344 //____________________________________________________________________________..
0345 void ThirdJetSpectra::Print(const std::string &what) const
0346 {
0347     std::cout << "ThirdJetSpectra::Print(const std::string &what) const Printing info for " << what << std::endl;
0348     //I will need to normalize the decorreclation average after the fact, but thats ok 
0349     std::cout<<"The number of events with at least three final state jets of pt> 3 GeV " <<n_three <<" and the number of events with an ID'ed dijet pair " <<n_evts <<std::endl;
0350     std::cout<<"looked at " <<n_try<<" events" <<std::endl; 
0351     TFile* f1=new TFile(Form("Third_Jet_Contrib_%s.root", segment_numb.c_str()), "RECREATE");
0352     first_jet_pt->Write();
0353     first_jet_E->Write();
0354     first_jet_eta->Write();
0355     first_jet_phi->Write(); 
0356     second_jet_pt->Write();
0357     second_jet_E->Write();
0358     second_jet_eta->Write();
0359     second_jet_phi->Write();    
0360     third_jet_pt->Write();
0361     third_jet_E->Write();
0362     third_jet_eta->Write();
0363     third_jet_phi->Write(); 
0364 
0365     dphi_12->Write();
0366     dphi_13->Write();
0367     dphi_23->Write();
0368     dphi_123->Write();
0369     pt_123->Write();
0370 
0371     decorr->Write();
0372     decorr_strict->Write();
0373     third_jet_pt_dec->Write();
0374     third_jet_pt_dec_strict->Write();
0375     third_jet_pt_dec_n->Write();
0376     third_jet_pt_dec_strict_n->Write();
0377 
0378     xj->Write();
0379     xj_strict->Write();
0380     xj_onl->Write();
0381     xj_strict_onl->Write();
0382     xj_12->Write();
0383     xj_13->Write();
0384     xj_23->Write();
0385     xj_strict_12->Write();
0386     xj_strict_13->Write();
0387     xj_strict_23->Write();
0388 
0389     f1->Write();
0390     f1->Close();    
0391 }
[ 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!