LargeRLCaloENC/src/DijetEventCuts.h

0001 #ifndef __DIJETEVENTCUTS_H__
0002 #define __DIJETEVENTCUTS_H__
0003 //Fun4All Related
0004 #include <fun4all/Fun4AllBase.h>
0005 #include <jetbase/JetContainer.h>
0006 #include <jetbase/JetContainerv1.h>
0007 #include <jetbase/JetMapv1.h>
0008 #include <jetbase/JetMap.h>
0009 #include <jetbase/Jetv1.h>
0010 #include <jetbase/Jetv2.h>
0011 #include <jetbase/Jet.h>
0012
0013 #include <vector>
0014 #include <string>
0015 #include <math.h>
0016 #include <array>
0017 #include <TTree.h>
0018 #define PI 3.14159265358979323464
0019
0020 class Jet;
0021 class JetContainer;
0022
0023 class DijetEventCuts{
0024     //maybe this gets some histos added for easier QA-ing of cut safety
0025     public:
0026         DijetEventCuts(float lpt=12., float slpt=7., float det=0.7, float dph=2.94,float maxpct=0.9, bool dj=true, bool ne=false, std::string radius="r04" ):
0027                 leadingpt(lpt),
0028                 subleadingpt(slpt),
0029                 etaedge(det), //keep full jet in calo
0030                 deltaphi(dph), //two hcal towers
0031                 maxOHCAL(maxpct),
0032                 isdijet(dj),
0033                 negativeEnergy(ne)
0034         {
0035             JetCuts=new TTree(Form("cuts_%s", radius.c_str()) , Form("Jet Cut Kinematics for QA on the Jet Event for jet with radius %s", radius.c_str()));
0036             JetCuts->Branch("N_Jets", &m_nJets, "Number of Jets in Container/I");
0037             JetCuts->Branch("z_vtx", &m_zvtx, "z vertex position/F");
0038             JetCuts->Branch("Lead_pt", &m_lpt, "p_{T} of Leading Jet/F");
0039             JetCuts->Branch("Subleading_pt", &m_slpt, "p_{T} of subleading jet (pair if dijet pair exists)/F");
0040             JetCuts->Branch("Lead_eta", &m_etal, "#eta of leading jet center/F");
0041             JetCuts->Branch("Subleading_eta", &m_etasl, "#eta subleading jet center/F");
0042             JetCuts->Branch("delta_phi", &m_deltaphi, "m_deltaphi/F");
0043             JetCuts->Branch("ohcal_ratio", &m_ohcalrat, "m_ohcalrat/F");
0044
0045             JetCuts->Branch("lead_ohcal_ratio", &m_leadER, "m_leadER/F");
0046             JetCuts->Branch("subleading_ohcal_ratio", &m_subleadingER, "m_leadER/F");
0047             JetCuts->Branch("lead_emcal_ratio", &m_leadER_E, "m_leadER/F");
0048             JetCuts->Branch("subleading_emcal_ratio", &m_subleadingER_E, "m_leadER/F");
0049             JetCuts->Branch("isDijet", &m_isdijet, "m_isDijet/O");
0050             JetCuts->Branch("negE", &m_hasnege, "m_hasnege/O");
0051             JetCuts->Branch("Il_E", &m_ile, "Moment of Inertia of leading jet/F");
0052             JetCuts->Branch("Isl_E", &m_isle, "Moment of Inertia of subleading jet/F");
0053             JetCuts->Branch("passes", &passesCut, "passesCut/O");
0054         };
0055         TTree* JetCuts; //This is a QA testing ttree to allow for immediate QA
0056         bool passesTheCut(JetContainer* eventjets, std::vector<float> ohcal_ratio_jets, std::vector<float> emcal_ratio_jets, float hcalratio, std::array<float,3> vertex, std::vector<float> i_e){
0057
0058             //check if the particular event passed the cut
0059             bool good=true;
0060             m_ohcalrat=hcalratio;
0061             if(hcalratio < 0 ) good=false;
0062             if(hcalratio >0.95 /* maxOHCAL*/) good=false;
0063             if(abs(vertex[2]) > 30 ) good=false; //cut on z=30 vertex
0064             m_zvtx=vertex[2];
0065             float leadjetpt=0., subleadjetpt=0.;
0066
0067             float leadingenergyratio=0., subleadingenergyratio=0.; //ohcal energy ratio
0068             float leadingenergyratio_E=0., subleadingenergyratio_E=0.; //emcal energy ratio
0069             bool haspartner=false;
0070                 Jet* leadjet=NULL, *subleadjet=NULL;
0071             int index=0; //associate the jet with its energy ratios
0072
0073                 for(auto j: *eventjets){
0074                 float pt=j->get_pt();
0075                 if(pt > leadjetpt){
0076                     if(leadjet)subleadjet=leadjet;
0077                     leadjet=j;
0078                         subleadjetpt=leadjetpt;
0079                         leadjetpt=pt;
0080
0081                     leadingenergyratio=ohcal_ratio_jets.at(index);
0082                     leadingenergyratio_E=emcal_ratio_jets.at(index);
0083                     m_ile=i_e.at(index);
0084
0085                        }
0086                 if(!negativeEnergy){
0087                     if(j->get_e() < 0){
0088                         m_hasnege=true;
0089                             good=false;
0090                     }
0091                     }
0092
0093                 index++;
0094
0095                 }
0096             if(leadjetpt < leadingpt) good=false;
0097             if(!leadjet) good=false;
0098             else{
0099             leadeta=leadjet->get_eta();
0100             m_etal=leadeta;
0101
0102             m_leadER=leadingenergyratio;
0103             m_leadER_E=leadingenergyratio_E;
0104                 leadphi=leadjet->get_phi();
0105             if(abs(leadeta) > etaedge ) good=false; //getting rid of events that have the leading jet outside of acceptance region
0106             if( m_leadER > maxOHCAL ) good=false;
0107             for(auto j: *eventjets){
0108                 float phi=j->get_phi();
0109                 int index=0;
0110
0111                 if(abs(phi-leadphi) > deltaphi && abs(phi-leadphi) <= PI+0.2){
0112                         subleadjet=j;
0113                     subleadjetpt=j->get_pt();
0114                     haspartner=true;
0115
0116                     subleadingenergyratio=ohcal_ratio_jets.at(index);
0117                     subleadingenergyratio_E=emcal_ratio_jets.at(index);
0118                     m_isle=i_e.at(index);
0119
0120                 break;}
0121                 index++;
0122
0123             }
0124             if(subleadjet){
0125                 float sldeta=subleadjet->get_eta();
0126                 if(abs(sldeta) > etaedge) good=false;
0127                 m_etasl=sldeta;
0128                 m_deltaphi=subleadjet->get_phi()-leadjet->get_phi();
0129
0130                 m_subleadingER=subleadingenergyratio;
0131                 m_subleadingER_E=subleadingenergyratio_E;
0132                 if(m_subleadingER > maxOHCAL) good=false;
0133
0134             }
0135
0136 alse;
0137             }
0138             //if(subleadjetpt < subleadingpt || !haspartner) good=false;
0139             passesCut=good;
0140             m_isdijet=haspartner;
0141             m_nJets=eventjets->size();
0142             m_lpt=leadjetpt;
0143             m_slpt=subleadjetpt;
0144             JetCuts->Fill();
0145             return passesCut;
0146         }
0147         float getLeadPhi(){ return leadphi;}
0148         float getLeadEta(){ return leadeta;}
0149         void getDijets(JetContainer* event_jets, std::vector<std::array<float, 3>>* dijet_sets)
0150         {
0151             std::vector<std::pair<Jet*, Jet*>> dijet_pairs;
0152             for(auto jet1=event_jets->begin(); jet1 != event_jets->end(); ++jet1){
0153                 float eta1=(*jet1)->get_eta(), phi1=(*jet1)->get_phi();
0154                 if(abs(eta1) > etaedge) continue; //keep R=0.4 jets in sPHENIX acceptance
0155                 if((*jet1)->get_pt() < 1.) continue; //reject jets below 1 GeV
0156                 for(auto jet2=jet1; jet2 != event_jets->end(); ++jet2){
0157                     if((*jet1) == (*jet2)) continue; //just double check
0158                     if((*jet2)->get_pt() < etaedge) continue;
0159                     float eta2=(*jet2)->get_eta(), phi2=(*jet2)->get_phi();
0160                     if(abs(eta2) > 0.7)continue;
0161                     if(abs(phi1 - phi2) >= deltaphi)dijet_pairs.push_back(std::make_pair(*jet1, *jet2));
0162                 }
0163             }
0164             for(auto dj:dijet_pairs){
0165                 float pt1=dj.first->get_pt(), pt2=dj.second->get_pt();
0166                 if(pt1 < pt2){
0167                     float ptt=pt1;
0168                     pt1=pt2;
0169                     pt2=ptt;
0170                 }
0171                 float xj=pt2/pt1, Ajj=(pt1-pt2)/(pt1+pt2);
0172                 std::array<float, 3> dijet_kin={pt1, Ajj, xj};
0173                 dijet_sets->push_back(dijet_kin);
0174             }
0175             return;
0176         }
0177     private:
0178
0179         float leadingpt=0.;
0180         float subleadingpt=0.;
0181         float etaedge=0.;
0182         float deltaphi=0.;
0183         float maxOHCAL=0.;
0184         bool isdijet=false;
0185         bool negativeEnergy=false;
0186
0187         bool passesCut=false;
0188         int m_nJets=0;
0189         float m_zvtx=0.;
0190         float m_lpt=0. ;
0191         float m_slpt=0.;
0192         float m_etal=0.;
0193         float m_etasl=0.;
0194
0195         float m_ile=0.;
0196         float m_isle=0.;
0197         float m_deltaphi=0.;
0198         float m_ohcalrat=0.;
0199         float leadphi=0.;
0200         float leadeta=0.;
0201
0202         float m_subleadingER=0.;
0203         float m_leadER=0.;
0204         float m_subleadingER_E=0.;
0205         float m_leadER_E=0.;
0206         bool m_isdijet=false;
0207         bool m_hasnege=false;
0208 };
0209
0210 #endif