Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-05-15 07:41:25 

0001 //==============================================================================
0002 
0003 //==============================================================================
0004 
0005 #include "Pythia8/Pythia.h"
0006 #include "TTree.h"
0007 #include "TFile.h"
0008 #include "TMath.h"
0009 
0010 #include <vector>
0011 #include <utility>
0012 #include <map>
0013 
0014 #include "fastjet/ClusterSequence.hh"
0015 #include <iostream>
0016 
0017 #include "Pythia8Plugins/HepMC3.h"
0018 
0019 #include "PythiaEvent.h"
0020 #include "HistogramsPythia.h"
0021 
0022 #define PR(x) std::cout << #x << " = " << (x) << std::endl;
0023 
0024 using namespace fastjet;
0025 using namespace Pythia8; 
0026 using namespace std;
0027 
0028 double costhetastar(int, int, const Event&);
0029 bool isInAcceptance(int, const Event&);  // acceptance filter
0030 int MakeEvent(Pythia *pythia, PythiaEvent *eventStore, int iev, bool writeTree = false);            // event handler (analyze event and
0031                                          // stores data in tuple)
0032 int findFinalElectron(const Event& event);
0033 
0034 void FindPartonJet(vector<fastjet::PseudoJet> jets, Particle parton1, Particle parton2, int &jetid1, int &jetid2);
0035 double dR(fastjet::PseudoJet jet, Particle parton);
0036 
0037 int isInHcals(Particle part);
0038 int isInHcals(double eta);
0039 bool isInTracker(double eta);
0040 //int isInTracker(Particle part);
0041 
0042 int FillHCals(TH2F *hist, TH1F *hEne, TH1F *hEneDenom, bool &anyHcal_jets);
0043 int FillHCalsJets(TH2F *hist, vector<fastjet::PseudoJet> jets);
0044 int FillHCalsJetsShare(TH2F *hist, vector<fastjet::PseudoJet> jets);
0045 
0046 void SortPairs(std::vector<pair<int, double>> &input, std::vector<pair<int, double>> &output);
0047 
0048 bool compFunc(pair<int, double> a, pair<int, double> b);
0049 void PrintVec(std::vector<pair<int, double>> input);
0050 
0051 int main(int argc, char* argv[]) {
0052     
0053     if (argc != 3) {
0054         cout << "Usage: " << argv[0] << " runcard  outfile" << endl;
0055         return 2;
0056     }
0057     char* runcard  = argv[1];
0058     char* outfile = argv[2];
0059     //const char* xmlDB    = "/users/PAS2524/lkosarz/Pythia/pythia8312/share/Pythia8/xmldoc";
0060     const char* xmlDB    = "/opt/local/share/Pythia8/xmldoc";
0061     
0062     bool WriteHepMC = true;
0063     bool WriteTree = false;
0064 
0065     // https://pythia.org//latest-manual/examples/main345.html
0066 
0067     //
0068     //  Create instance of Pythia 
0069     //
0070     //Pythia pythia(xmlDB); // the default parameters are read from xml files
0071                           // stored in the xmldoc directory. This includes
0072                           // particle data and decay definitions.
0073     
0074     Pythia *pythia = new Pythia("/opt/local/share/Pythia8/xmldoc");
0075 
0076 
0077     // Shorthand for (static) settings
0078     Settings& settings = pythia->settings;
0079     
0080     //  Read in runcard (should be star_hf_tune_v1.0.cmd)
0081     pythia->readFile(runcard);
0082     cout << "Runcard '" << runcard << "' loaded." << endl;
0083     
0084 
0085     TFile *treeFile  = new TFile(Form("%s_tree.root", outfile),"RECREATE");
0086     TFile *histFile  = new TFile(Form("%s_hist.root", outfile),"RECREATE");
0087 
0088     treeFile->cd();
0089 
0090     TTree *eventTree = new TTree("eventTree", "Event tree");
0091     PythiaEvent *eventStore = new PythiaEvent();
0092     eventTree->Branch("eventTree", &eventStore, 32000, 1);
0093 
0094 
0095 
0096 
0097     //
0098     //  Retrieve number of events and other parameters from the runcard.
0099     //  We need to deal with those settings ourself. Getting
0100     //  them through the runcard just avoids recompiling.
0101     //
0102     long  maxNumberOfEvents = settings.mode("Main:numberOfEvents");
0103     //int  nList     = settings.mode("Main:numberToList");
0104     int  nList = 10;
0105     //int  nShow     = settings.mode("Main:timesToShow");
0106     int  nShow = 10;
0107     int  maxErrors = settings.mode("Main:timesAllowErrors");
0108     //bool showCS    = settings.flag("Main:showChangedSettings");
0109     //bool showAS    = settings.flag("Main:showAllSettings");
0110     int  pace = maxNumberOfEvents/nShow;
0111  
0112     
0113     //  Initialize Pythia, ready to go
0114     pythia->init();
0115     HepMC3::WriterAscii hepmcWriter(Form("%s.hepmc3", outfile));
0116     HepMC3::Pythia8ToHepMC3 toHepMC;
0117     
0118     // List changed or all data
0119     //if (showCS) settings.listChanged();
0120     //if (showAS) settings.listAll();
0121     settings.listChanged();
0122     settings.listAll();
0123     
0124 
0125     
0126     //IO_GenEvent ascii_io(Form("%s.hepmc3", outfile));
0127     //HepMC::WriterRootTree WriterRootfile("file.root")
0128     //GenEvent hepmcevt;
0129 
0130 
0131     histFile->cd();
0132     CreateHistograms();
0133 
0134     //--------------------------------------------------------------
0135     //  Event loop
0136     //--------------------------------------------------------------
0137     int ievent = 0;
0138     int iErrors = 0;
0139     
0140     while (ievent < maxNumberOfEvents) {
0141         
0142         if (!pythia->next()) {
0143             if (++iErrors < maxErrors) continue;
0144             cout << "Error: too many errors in event generation - check your settings & code" << endl;
0145             break;
0146         }
0147 
0148         h_Events_cuts->Fill(0);
0149 
0150         h_XsecGen->Fill(pythia->info.sigmaGen());
0151         h_XsecGen_err->Fill(pythia->info.sigmaGen(), pythia->info.sigmaErr());
0152 
0153         //--------------------
0154 
0155         bool isDiffractive = pythia->info.isDiffractiveA() || pythia->info.isDiffractiveB();
0156         bool isHardDiffractive = pythia->info.isHardDiffractiveA() || pythia->info.isHardDiffractiveB();
0157 
0158         //if (isDiffractive) cout<<"Diffractive"<<endl;
0159         //if (isHardDiffractive) cout<<"Hard diffractive"<<endl;
0160 
0161         if (isDiffractive) h_Events_cuts->Fill(1);
0162         if (isHardDiffractive) h_Events_cuts->Fill(2);
0163 
0164         //if(!(isDiffractive || isHardDiffractive)) continue;
0165 
0166         int phMode = pythia->info.photonMode();
0167 
0168         //-------------------------
0169 
0170         h_XsecSel->Fill(pythia->info.sigmaGen());
0171         h_XsecSel_err->Fill(pythia->info.sigmaGen(), pythia->info.sigmaErr());
0172 
0173         MakeEvent(pythia, eventStore, ievent, WriteTree);  // in MakeEvent we deal with the whole event and return
0174 
0175         if(WriteTree)
0176         {
0177             eventTree->Fill();
0178             eventStore->Clear();
0179         }
0180         ievent++;
0181 
0182         if (ievent%pace == 0) {
0183             cout << "# of events generated = " << ievent << endl;
0184         }
0185         
0186         // List first few events.
0187         if (ievent < nList) {
0188             pythia->info.list();
0189             pythia->process.list();
0190             pythia->event.list(false, true, 3);
0191         }
0192 
0193 
0194         if (WriteHepMC) {
0195             HepMC3::GenEvent evt;
0196             toHepMC.fill_next_event(*pythia, &evt);
0197             hepmcWriter.write_event(evt);
0198         }
0199 
0200         //GenEvent* hepmcevt = new HepMC::GenEvent();
0201         //toHepMC.fill_next_event( *pythia.Pythia8(), hepmcevt);
0202         //ascii_io << hepmcevt;
0203         //delete hepmcevt;
0204 
0205     }
0206     
0207     //--------------------------------------------------------------
0208     //  Finish up
0209     //--------------------------------------------------------------
0210     //pythia.statistics();
0211     pythia->stat();
0212 
0213     cout << "Writing files" << endl;
0214 
0215     histFile->cd();
0216     histFile->Write();
0217 
0218     treeFile->cd();
0219     if(WriteTree) eventTree->Write();
0220     //treeFile->Write();
0221 
0222     histFile->Close();
0223     treeFile->Close();
0224     hepmcWriter.close();
0225 
0226     cout << "Finish!" << endl;
0227     
0228     return 0;
0229 }
0230 
0231 //
0232 //  Event analysis
0233 //
0234 int MakeEvent(Pythia *pythia, PythiaEvent *eventStore, int iev, bool writeTree)
0235 {
0236     Event &event = pythia->event;
0237 
0238     //cout<<"NEW EVENT ---------------------------"<<endl;
0239 
0240     //std::cout<<"simulating event "<<iev<<std::endl;
0241 
0242 
0243     if(pythia->info.isDiffractiveA()) h_Events_Diffractive->Fill(0);
0244     if(pythia->info.isDiffractiveB()) h_Events_Diffractive->Fill(1);
0245     if(pythia->info.isHardDiffractiveA()) h_Events_Diffractive->Fill(2);
0246     if(pythia->info.isHardDiffractiveB()) h_Events_Diffractive->Fill(3);
0247 
0248     
0249     hist_eta_energy_tmp->Reset();
0250     hist_eta_energy_denom_tmp->Reset();
0251 
0252     //if (event[1].id() == 11)
0253     //int eleid = findFinalElectron(event);
0254 
0255     // create a jet definition:
0256     // a jet algorithm with a given radius parameter
0257     //----------------------------------------------------------
0258     double R = 1.0;
0259     double p = -1.0;
0260     fastjet::JetDefinition jet_def(fastjet::antikt_algorithm, R);
0261     //fastjet::JetDefinition jet_def(fastjet::ee_kt_algorithm);
0262     //fastjet::JetDefinition jet_def(fastjet::ee_genkt_algorithm, R, p);
0263     //fastjet::JetDefinition jet_def_meas(fastjet::ee_kt_algorithm);
0264 
0265     vector<fastjet::PseudoJet> input_particles;
0266     vector<fastjet::PseudoJet> input_particles_meas;
0267     vector<fastjet::PseudoJet> input_particles_meas_no_nHCal;
0268 
0269     float EsumF = 0.0;
0270     double EsumD = 0.0;
0271 
0272     bool has_nHcalActivity = false;
0273 
0274     int nPartFinal = 0;
0275     int nPartonsOut = 0;
0276     int iParton_1 = 0;
0277     int iParton_2 = 0;
0278 
0279     bool Parton_1_inAcc = false;
0280     bool Parton_2_inAcc = false;
0281 
0282     int nJetsInAcc = 0;
0283 
0284     int nPion_p = 0;
0285     int nPion_n = 0;
0286     int nKaon_p = 0;
0287     int nKaon_n = 0;
0288     int nProton_p = 0;
0289     int nProton_n = 0;
0290     int nElectron_p = 0;
0291     int nElectron_n = 0;
0292 
0293     int nNeutron = 0;
0294     int nGamma = 0;
0295 
0296 
0297     // Four-momenta of proton, electron, virtual photon/Z^0/W^+-.
0298     Vec4 pProton = event[1].p();
0299     Vec4 peIn    = event[2].p();
0300     Vec4 pPhoton = event[4].p();
0301 
0302     //Vec4 peIn    = event[4].p();
0303     //Vec4 peOut   = event[6].p();
0304     //Vec4 pPhoton = peIn - peOut;
0305 
0306     // Q2, W2, Bjorken x, y.
0307     double Q2    = - pPhoton.m2Calc();
0308     double W2    = (pProton + pPhoton).m2Calc();
0309     double x     = Q2 / (2. * pProton * pPhoton);
0310     double y     = (pProton * pPhoton) / (pProton * peIn);
0311 
0312 
0313     for (int i = 0; i < event.size(); i++) {
0314 
0315         Particle part = event[i];
0316 
0317          // count outgoing partons
0318         if(part.status()==-23 || part.status()==-24)
0319         {
0320             nPartonsOut++;
0321 
0322             h_Partons_status->Fill(part.status());
0323 
0324             h_Parton_eta_p->Fill(part.eta(), part.pAbs());
0325             h_Parton_eta_pT->Fill(part.eta(), part.pT());
0326             h_Parton_eta_E->Fill(part.eta(), part.e());
0327 
0328             if(nPartonsOut == 1) iParton_1 = i;
0329             if(nPartonsOut == 2) iParton_2 = i;
0330 
0331             if(part.id()>0) h_Partons_types->Fill(part.id());
0332             if(part.id()<0) h_Partons_types_anti->Fill(abs(part.id()));
0333 
0334             if((part.eta()>=-4.14 && part.eta()<=4.2) && nPartonsOut == 1) Parton_1_inAcc = true;
0335             if((part.eta()>=-4.14 && part.eta()<=4.2) && nPartonsOut == 2) Parton_2_inAcc = true;
0336         }
0337 
0338         if(nPartonsOut == 2)
0339         {
0340             h_Partons_eta->Fill(event[iParton_1].eta(), event[iParton_2].eta());
0341             h_Partons_phi->Fill(event[iParton_1].phi(), event[iParton_2].phi());
0342             h_Partons_p->Fill(event[iParton_1].pAbs(), event[iParton_2].pAbs());
0343             h_Partons_pT->Fill(event[iParton_1].pT(), event[iParton_2].pT());
0344 
0345             h_Parton_x_eta->Fill(x, event[iParton_1].eta());
0346             h_Parton_x_eta->Fill(x, event[iParton_2].eta());
0347 
0348             h_Parton_y_eta->Fill(y, event[iParton_1].eta());
0349             h_Parton_y_eta->Fill(y, event[iParton_2].eta());
0350 
0351             h_Parton_x_eta1->Fill(x, event[iParton_1].eta());
0352             h_Parton_x_eta2->Fill(x, event[iParton_2].eta());
0353 
0354             h_Parton_y_eta1->Fill(y, event[iParton_1].eta());
0355             h_Parton_y_eta2->Fill(y, event[iParton_2].eta());
0356         }
0357 
0358         // accept final particles only
0359         if(!part.isFinal()) continue;
0360         // ePIC acceptance only
0361         if(part.eta()<-4.14 ||  part.eta()>4.2) continue;
0362         if(!(11>part.status() || part.status()>19)) continue; // exclude beam particles
0363 
0364         nPartFinal++;
0365 
0366 
0367         if(part.id() == 211) nPion_p++;
0368         if(part.id() == -211) nPion_n++;
0369         if(part.id() == 321) nKaon_p++;
0370         if(part.id() == -321) nKaon_n++;
0371         if(part.id() == 2212) nProton_p++;
0372         if(part.id() == -2212) nProton_n++;
0373         if(part.id() == -11) nElectron_p++;
0374         if(part.id() == 11) nElectron_n++;
0375 
0376         if(part.id() == 2112) nNeutron++;
0377         if(part.id() == 22) nGamma++;
0378 
0379         h_Particle_eta->Fill(part.eta());
0380         h_Particle_eta_wE->Fill(part.eta(), part.e());
0381 
0382         h_Particle_eta_p->Fill(part.eta(), part.pAbs());
0383         h_Particle_eta_pT->Fill(part.eta(), part.pT());
0384         h_Particle_eta_E->Fill(part.eta(), part.e());
0385 
0386 
0387         // eta, momentum
0388         if(part.id() == 211) h_Particle_pion_p_eta_p->Fill(part.eta(), part.pAbs());
0389         if(part.id() == -211) h_Particle_pion_n_eta_p->Fill(part.eta(), part.pAbs());
0390         if(part.id() == 321) h_Particle_Kaon_p_eta_p->Fill(part.eta(), part.pAbs());
0391         if(part.id() == -321) h_Particle_Kaon_n_eta_p->Fill(part.eta(), part.pAbs());
0392         if(part.id() == 2212) h_Particle_proton_p_eta_p->Fill(part.eta(), part.pAbs());
0393         if(part.id() == -2212) h_Particle_proton_n_eta_p->Fill(part.eta(), part.pAbs());
0394         if(part.id() == -11) h_Particle_Electron_p_eta_p->Fill(part.eta(), part.pAbs());
0395         if(part.id() == 11) h_Particle_Electron_n_eta_p->Fill(part.eta(), part.pAbs());
0396 
0397         if(part.id() == 2112) h_Particle_Neutron_eta_p->Fill(part.eta(), part.pAbs());
0398         if(part.id() == 22) h_Particle_Gamma_eta_p->Fill(part.eta(), part.pAbs());
0399 
0400         // eta, transverse momentum pT
0401         if(part.id() == 211) h_Particle_pion_p_eta_pT->Fill(part.eta(), part.pT());
0402         if(part.id() == -211) h_Particle_pion_n_eta_pT->Fill(part.eta(), part.pT());
0403         if(part.id() == 321) h_Particle_Kaon_p_eta_pT->Fill(part.eta(), part.pT());
0404         if(part.id() == -321) h_Particle_Kaon_n_eta_pT->Fill(part.eta(), part.pT());
0405         if(part.id() == 2212) h_Particle_proton_p_eta_pT->Fill(part.eta(), part.pT());
0406         if(part.id() == -2212) h_Particle_proton_n_eta_pT->Fill(part.eta(), part.pT());
0407         if(part.id() == -11) h_Particle_Electron_p_eta_pT->Fill(part.eta(), part.pT());
0408         if(part.id() == 11) h_Particle_Electron_n_eta_pT->Fill(part.eta(), part.pT());
0409 
0410         if(part.id() == 2112) h_Particle_Neutron_eta_pT->Fill(part.eta(), part.pT());
0411         if(part.id() == 22) h_Particle_Gamma_eta_pT->Fill(part.eta(), part.pT());
0412 
0413         // eta, energy
0414         if(part.id() == 211) h_Particle_Pion_p_eta_E->Fill(part.eta(), part.e());
0415         if(part.id() == -211) h_Particle_Pion_n_eta_E->Fill(part.eta(), part.e());
0416         if(part.id() == 321) h_Particle_Kaon_p_eta_E->Fill(part.eta(), part.e());
0417         if(part.id() == -321) h_Particle_Kaon_n_eta_E->Fill(part.eta(), part.e());
0418         if(part.id() == 2212) h_Particle_Proton_p_eta_E->Fill(part.eta(), part.e());
0419         if(part.id() == -2212) h_Particle_Proton_n_eta_E->Fill(part.eta(), part.e());
0420         if(part.id() == -11) h_Particle_Electron_p_eta_E->Fill(part.eta(), part.e());
0421         if(part.id() == 11) h_Particle_Electron_n_eta_E->Fill(part.eta(), part.e());
0422 
0423         if(part.id() == 2112) h_Particle_Neutron_eta_E->Fill(part.eta(), part.e());
0424         if(part.id() == 22) h_Particle_Gamma_eta_E->Fill(part.eta(), part.e());
0425 
0426 
0427         // read in input particles
0428         //----------------------------------------------------------
0429         input_particles.push_back(fastjet::PseudoJet(part.px(),part.py(),part.pz(),part.e()));
0430         //cout<<"input_particles add = "<<i<<endl;
0431 
0432         if(isInHcals(part.eta())>0)
0433         {
0434             if(part.isCharged() && isInTracker(part.eta()))
0435             {
0436                 input_particles_meas.push_back(fastjet::PseudoJet(part.px(),part.py(),part.pz(),part.e()));
0437                 //cout<<"input_particles_meas add charged = "<<i<<endl;
0438             }
0439             if(part.isNeutral())
0440             {
0441                 input_particles_meas.push_back(fastjet::PseudoJet(part.px(),part.py(),part.pz(),part.e()));
0442                 //cout<<"input_particles_meas add neutral = "<<i<<endl;
0443             }
0444 
0445             bool isIn_nHCal = -4.14<part.eta() && part.eta()<-1.18;
0446 
0447             if(part.isCharged() && isInTracker(part.eta()))
0448             {
0449                 input_particles_meas_no_nHCal.push_back(fastjet::PseudoJet(part.px(),part.py(),part.pz(),part.e()));
0450                 //cout<<"input_particles_meas_no_nHCal add charged = "<<i<<endl;
0451             }
0452             if(part.isNeutral() && !isIn_nHCal)
0453             {
0454                 input_particles_meas_no_nHCal.push_back(fastjet::PseudoJet(part.px(),part.py(),part.pz(),part.e()));
0455                 //cout<<"input_particles_meas_no_nHCal add neutral = "<<i<<endl;
0456             }
0457 
0458         }
0459 
0460 
0461         if(11>part.status() || part.status()>19) // exclude beam particles
0462         {
0463             hist_eta_energy_tmp->Fill(part.eta(), part.e());
0464 
0465             //cout<<"part E ="<<part.e()<<endl;
0466 
0467             EsumF += part.e();
0468             EsumD += part.e();
0469 
0470             //cout<<"EsumF ="<<EsumF<<endl;
0471             //cout<<"EsumD ="<<EsumD<<endl;
0472 
0473             for (int i = 1; i <= hist_eta_energy_denom_tmp->GetNbinsX(); ++i) {
0474 
0475                 double xbin_cent = hist_eta_energy_denom_tmp->GetXaxis()->GetBinCenter(i);
0476 
0477                 hist_eta_energy_denom_tmp->Fill(xbin_cent, part.e());
0478 
0479             }
0480         }
0481 
0482 
0483         // Store particle
0484         if(writeTree) eventStore->particles.push_back(part);
0485 
0486 
0487         // select nHCal acceptance
0488         //if(part.eta()<-4.03 ||  part.eta()>-1.27) continue;
0489         if(part.eta()<-4.14 ||  part.eta()>-1.18)   continue;
0490 
0491         has_nHcalActivity = true;
0492 
0493 
0494 
0495         /*
0496             //
0497             //  Get daughters
0498             //
0499             vector<int> daughters = event.daughterList(i);
0500             
0501             if (daughters.size() != 2) continue;
0502             int ielectron = daughters[0];
0503             int ipositron = daughters[1];
0504             
0505             if (abs(event[ielectron].id()) != 11) continue;
0506             if (abs(event[ipositron].id()) != 11) continue;
0507             
0508             if (event[ielectron].id() == -11) {
0509                 int k = ielectron;
0510                 ielectron = ipositron;
0511                 ipositron = k;
0512             }*/
0513             
0514            // if (!(isInAcceptance(ielectron, event) && isInAcceptance(ipositron, event))) continue;
0515             
0516             
0517     } // particles loop
0518 
0519 
0520     // run the jet clustering with the above jet definition
0521     //----------------------------------------------------------
0522     fastjet::ClusterSequence clust_seq(input_particles, jet_def);
0523     fastjet::ClusterSequence clust_seq_meas(input_particles_meas, jet_def);
0524     fastjet::ClusterSequence clust_seq_meas_no_nHCal(input_particles_meas_no_nHCal, jet_def);
0525 
0526 
0527     // get the resulting jets ordered in pt
0528     //----------------------------------------------------------
0529     double ptmin = 4.0; // 5.0 [GeV/c]
0530     double Emin = 0.0; // 5.0 [GeV/c]
0531     vector<fastjet::PseudoJet> inclusive_jets_unsorted;
0532     vector<fastjet::PseudoJet> measured_jets_unsorted;
0533     vector<fastjet::PseudoJet> measured_jets_no_nHCal_unsorted;
0534 
0535     vector<fastjet::PseudoJet> inclusive_jets;
0536     vector<fastjet::PseudoJet> measured_jets;
0537     vector<fastjet::PseudoJet> measured_jets_no_nHCal;
0538 
0539     //vector<fastjet::PseudoJet> inclusive_jets = sorted_by_pt(clust_seq.inclusive_jets(ptmin));
0540     //inclusive_jets = sorted_by_E(clust_seq.inclusive_jets(ptmin));
0541     //measured_jets = sorted_by_E(clust_seq_meas.inclusive_jets(ptmin));
0542 
0543     vector<fastjet::PseudoJet> inclusive_jets_precut = sorted_by_E(clust_seq.inclusive_jets(ptmin));
0544     vector<fastjet::PseudoJet> measured_jets_precut = sorted_by_E(clust_seq_meas.inclusive_jets(ptmin));
0545     vector<fastjet::PseudoJet> measured_jets_no_nHCal_precut = sorted_by_E(clust_seq_meas_no_nHCal.inclusive_jets(ptmin));
0546 
0547 
0548     // apply energy cut
0549     for (int i = 0; i < inclusive_jets_precut.size(); ++i)
0550     {
0551         fastjet::PseudoJet jet = inclusive_jets_precut[i];
0552 
0553         if(jet.E()<Emin) continue;
0554 
0555         if(jet.eta()>=-4.14 &&  jet.eta()<=4.2) nJetsInAcc++;
0556 
0557         inclusive_jets_unsorted.push_back(jet);
0558     }
0559 
0560     for (int i = 0; i < measured_jets_precut.size(); ++i)
0561     {
0562         fastjet::PseudoJet jet = measured_jets_precut[i];
0563 
0564         if(jet.E()<Emin) continue;
0565 
0566         measured_jets_unsorted.push_back(jet);
0567     }
0568 
0569     for (int i = 0; i < measured_jets_no_nHCal_precut.size(); ++i)
0570     {
0571         fastjet::PseudoJet jet = measured_jets_no_nHCal_precut[i];
0572 
0573         if(jet.E()<Emin) continue;
0574 
0575         measured_jets_no_nHCal_unsorted.push_back(jet);
0576     }
0577 
0578     // sort jets vs. energy
0579 /*
0580     inclusive_jets = sorted_by_E(inclusive_jets_unsorted);
0581     measured_jets = sorted_by_E(measured_jets_unsorted);
0582     measured_jets_no_nHCal = sorted_by_E(measured_jets_no_nHCal_unsorted);
0583     */
0584 
0585     // sort jets vs. pT
0586     inclusive_jets = sorted_by_pt(inclusive_jets_unsorted);
0587     measured_jets = sorted_by_pt(measured_jets_unsorted);
0588     measured_jets_no_nHCal = sorted_by_pt(measured_jets_no_nHCal_unsorted);
0589 
0590     //cout<<"inclusive_jets size = "<<inclusive_jets.size()<<endl;
0591     //cout<<"measured_jets size = "<<measured_jets.size()<<endl;
0592 
0593 
0594     h_Events->Fill(0);
0595 
0596     h_Events_types->Fill(0);
0597     if((Parton_1_inAcc && !Parton_2_inAcc) || (Parton_2_inAcc && !Parton_1_inAcc)) h_Events_types->Fill(1);
0598     if(Parton_1_inAcc && Parton_2_inAcc) h_Events_types->Fill(2);
0599     if(nJetsInAcc == 1) h_Events_types->Fill(3);
0600     if(nJetsInAcc == 2) h_Events_types->Fill(4);
0601     if(nJetsInAcc >= 2) h_Events_types->Fill(5);
0602 
0603     if(Parton_1_inAcc && Parton_2_inAcc) h_Events_cuts->Fill(3);
0604 
0605     h_Events_nPartonsOut->Fill(nPartonsOut);
0606 
0607     // Store event
0608     if(writeTree)
0609     {
0610         eventStore->eventId = iev;
0611         eventStore->nParticlesFinal = nPartFinal;
0612         eventStore->Q2 = Q2;
0613         eventStore->x = x;
0614         eventStore->y = y;
0615     }
0616 
0617 
0618     h_Event_xQ2->Fill(x, Q2);
0619     h_Event_yQ2->Fill(y, Q2);
0620     h_Event_xy->Fill(x, y);
0621 
0622 
0623     if(has_nHcalActivity)
0624     {
0625         h_Event_nHCal_xQ2->Fill(x, Q2);
0626         h_Event_nHCal_yQ2->Fill(y, Q2);
0627         h_Event_nHCal_xy->Fill(x, y);
0628     }
0629 
0630     h_Event_nPart_final->Fill(nPartFinal);
0631 
0632     h_Event_nPion_p->Fill(nPion_p);
0633     h_Event_nPion_n->Fill(nPion_n);
0634     h_Event_nKaon_p->Fill(nKaon_p);
0635     h_Event_nKaon_n->Fill(nKaon_n);
0636     h_Event_nProton_p->Fill(nProton_p);
0637     h_Event_nProton_n->Fill(nProton_n);
0638     h_Event_nElectron_p->Fill(nElectron_p);
0639     h_Event_nElectron_n->Fill(nElectron_n);
0640 
0641     h_Event_nNeutron->Fill(nNeutron);
0642     h_Event_nGamma->Fill(nGamma);
0643 
0644     h_Event_nJets->Fill(inclusive_jets.size());
0645     h_Event_nJets_meas->Fill(measured_jets.size());
0646     h_Event_nJets_meas_no_nHCal->Fill(measured_jets_no_nHCal.size());
0647 
0648     // summary
0649     bool anyHcal_jets = false;
0650     bool anyHcal_jets_meas = false;
0651     bool anyHcal_jets_meas_no_nHCal = false;
0652     bool bHCalJet = false;
0653     bool bHCalJet_meas = false;
0654     int nHCal_jets = 0;
0655     int nHCal_jets_meas = 0;
0656 
0657     //int nHCal_jets = FillHCals(h_Event_HCal_jets, hist_eta_energy_tmp, hist_eta_energy_denom_tmp, anyHcal_jets);
0658     FillHCalsJets(h_Event_HCal_jets, inclusive_jets);
0659     FillHCalsJets(h_Event_HCal_jets_meas, measured_jets);
0660     FillHCalsJets(h_Event_HCal_jets_meas_no_nHCal, measured_jets_no_nHCal);
0661 
0662     FillHCalsJetsShare(h_Event_HCal_jets_meas_full, measured_jets);
0663 
0664 
0665     for (unsigned int i = 0; i < inclusive_jets.size(); i++) {
0666 
0667         fastjet::PseudoJet jet = inclusive_jets[i];
0668 
0669         int HCal_id = 0; // none = 0, nHCal = 1, bHCal = 2, LFHCAL = 3
0670 
0671         // nHCal
0672         if(-4.14 < jet.eta() && jet.eta() < -1.18)
0673         {
0674             anyHcal_jets = true;
0675             nHCal_jets++;
0676             HCal_id = 1;
0677         }
0678 
0679         // bHCal
0680         if(-1.1 < jet.eta() && jet.eta() < 1.1)
0681         {
0682             anyHcal_jets = true;
0683             bHCalJet = true;
0684             HCal_id = 2;
0685         }
0686 
0687         // LFHCal
0688         if(1.2 < jet.eta() && jet.eta() < 4.2)
0689         {
0690             anyHcal_jets = true;
0691             HCal_id = 3;
0692         }
0693 
0694 
0695         double jet_p = sqrt(jet.pt2()+jet.pz()*jet.pz());
0696 
0697         h_Jet_nPart->Fill(jet.constituents().size());
0698         h_Jet_mass->Fill(jet.m());
0699         //h_Jet_charge->Fill(jet.m());
0700         h_Jet_E->Fill(jet.E());
0701         h_Jet_p->Fill(jet_p);
0702         h_Jet_pT->Fill(jet.pt());
0703         h_Jet_eta->Fill(jet.eta());
0704 
0705 
0706         for (unsigned int j = i; j < inclusive_jets.size(); j++) {
0707 
0708             if(i==j) continue;
0709 
0710             fastjet::PseudoJet jet2 = inclusive_jets[j];
0711 
0712             h_Jet_deta->Fill(jet.eta()-jet2.eta());
0713         }
0714 
0715         vector<PseudoJet> constituents = jet.constituents();
0716 
0717         for (int j = 0; j < constituents.size(); ++j) {
0718 
0719             if(bHCalJet) h_Jet_bHCal_part_eta->Fill(constituents[j].eta());
0720             h_Jet_HCal_part_eta->Fill(constituents[j].eta(), HCal_id);
0721         }
0722 
0723     } // inclusive_jets
0724 
0725     for (unsigned int i = 0; i < measured_jets.size(); i++) {
0726 
0727         fastjet::PseudoJet jet = measured_jets[i];
0728 
0729         int HCal_meas_id = 0; // none = 0, nHCal = 1, bHCal = 2, LFHCAL = 3
0730 
0731         // nHCal
0732         if(-4.14 < jet.eta() && jet.eta() < -1.18)
0733         {
0734             anyHcal_jets_meas = true;
0735             nHCal_jets_meas++;
0736             HCal_meas_id = 1;
0737         }
0738 
0739         // bHCal
0740         if(-1.1 < jet.eta() && jet.eta() < 1.1)
0741         {
0742             anyHcal_jets_meas = true;
0743             bHCalJet_meas = true;
0744             HCal_meas_id = 2;
0745         }
0746 
0747         // LFHCal
0748         if(1.2 < jet.eta() && jet.eta() < 4.2)
0749         {
0750             anyHcal_jets_meas = true;
0751             HCal_meas_id = 3;
0752         }
0753 
0754 
0755         double jet_p = sqrt(jet.pt2()+jet.pz()*jet.pz());
0756 
0757         h_Jet_meas_nPart->Fill(jet.constituents().size());
0758         h_Jet_meas_mass->Fill(jet.m());
0759         //h_Jet_meas_charge->Fill(jet.m());
0760         h_Jet_meas_E->Fill(jet.E());
0761         h_Jet_meas_p->Fill(jet_p);
0762         h_Jet_meas_pT->Fill(jet.pt());
0763         h_Jet_meas_eta->Fill(jet.eta());
0764 
0765         for (unsigned int j = i; j < measured_jets.size(); j++) {
0766 
0767             if(i==j) continue;
0768 
0769             fastjet::PseudoJet jet2 = measured_jets[j];
0770 
0771             h_Jet_meas_deta->Fill(jet.eta()-jet2.eta());
0772         }
0773 
0774 
0775         vector<PseudoJet> measured_constituents = jet.constituents();
0776 
0777         for (int j = 0; j < measured_constituents.size(); ++j) {
0778 
0779             if(bHCalJet) h_Jet_meas_bHCal_part_eta->Fill(measured_constituents[j].eta());
0780             h_Jet_meas_HCal_part_eta->Fill(measured_constituents[j].eta(), HCal_meas_id);
0781         }
0782     } // meas_jets
0783 
0784 
0785 
0786     for (unsigned int i = 0; i < measured_jets_no_nHCal.size(); i++) {
0787 
0788         fastjet::PseudoJet jet = measured_jets_no_nHCal[i];
0789 
0790         //int HCal_meas_id = 0; // none = 0, nHCal = 1, bHCal = 2, LFHCAL = 3
0791 /*
0792         // nHCal
0793         if(-4.14 < jet.eta() && jet.eta() < -1.18)
0794         {
0795             anyHcal_jets_meas = true;
0796             nHCal_jets_meas++;
0797             HCal_meas_id = 1;
0798         }*/
0799 
0800         // bHCal
0801         if(-1.1 < jet.eta() && jet.eta() < 1.1)
0802         {
0803             anyHcal_jets_meas_no_nHCal = true;
0804             //bHCalJet_meas = true;
0805             //HCal_meas_id = 2;
0806         }
0807 
0808         // LFHCal
0809         if(1.2 < jet.eta() && jet.eta() < 4.2)
0810         {
0811             anyHcal_jets_meas_no_nHCal = true;
0812             //HCal_meas_id = 3;
0813         }
0814 
0815 
0816         double jet_p = sqrt(jet.pt2()+jet.pz()*jet.pz());
0817 
0818         h_Jet_meas_no_nHCal_nPart->Fill(jet.constituents().size());
0819         h_Jet_meas_no_nHCal_mass->Fill(jet.m());
0820         //h_Jet_meas_no_nHCal_charge->Fill(jet.m());
0821         h_Jet_meas_no_nHCal_E->Fill(jet.E());
0822         h_Jet_meas_no_nHCal_p->Fill(jet_p);
0823         h_Jet_meas_no_nHCal_pT->Fill(jet.pt());
0824         h_Jet_meas_no_nHCal_eta->Fill(jet.eta());
0825 
0826     } // measured_jets_no_nHCal
0827 
0828 
0829 
0830 
0831     // jets ------------------
0832 
0833     if(inclusive_jets.size() >= 2)
0834     {
0835 
0836         fastjet::PseudoJet jet1 = inclusive_jets[0];
0837         fastjet::PseudoJet jet2 = inclusive_jets[1];
0838 
0839         double jet1_p = sqrt(jet1.pt2()+jet1.pz()*jet1.pz());
0840         double jet2_p = sqrt(jet2.pt2()+jet2.pz()*jet2.pz());
0841 
0842         h_Jets_eta->Fill(jet1.eta(), jet2.eta());
0843         h_Jets_phi->Fill(jet1.phi(), jet2.phi());
0844         h_Jets_p->Fill(jet1_p, jet2_p);
0845         h_Jets_pT->Fill(jet1.pt(), jet2.pt());
0846         h_Jets_E->Fill(jet1.E(), jet2.E());
0847     }
0848 
0849     if(measured_jets.size() >= 2)
0850     {
0851         h_Jets_meas_eta->Fill(measured_jets[0].eta(), measured_jets[1].eta());
0852 
0853         fastjet::PseudoJet jet1 = measured_jets[0];
0854         fastjet::PseudoJet jet2 = measured_jets[1];
0855 
0856         double jet1_p = sqrt(jet1.pt2()+jet1.pz()*jet1.pz());
0857         double jet2_p = sqrt(jet2.pt2()+jet2.pz()*jet2.pz());
0858 
0859         h_Jets_meas_eta->Fill(jet1.eta(), jet2.eta());
0860         h_Jets_meas_phi->Fill(jet1.phi(), jet2.phi());
0861         h_Jets_meas_p->Fill(jet1_p, jet2_p);
0862         h_Jets_meas_pT->Fill(jet1.pt(), jet2.pt());
0863         h_Jets_meas_E->Fill(jet1.E(), jet2.E());
0864     }
0865 
0866     if(measured_jets_no_nHCal.size() >= 2)
0867     {
0868         h_Jets_meas_eta->Fill(measured_jets_no_nHCal[0].eta(), measured_jets_no_nHCal[1].eta());
0869 
0870         fastjet::PseudoJet jet1 = measured_jets_no_nHCal[0];
0871         fastjet::PseudoJet jet2 = measured_jets_no_nHCal[1];
0872 
0873         double jet1_p = sqrt(jet1.pt2()+jet1.pz()*jet1.pz());
0874         double jet2_p = sqrt(jet2.pt2()+jet2.pz()*jet2.pz());
0875 
0876         h_Jets_meas_no_nHCal_eta->Fill(jet1.eta(), jet2.eta());
0877         h_Jets_meas_no_nHCal_phi->Fill(jet1.phi(), jet2.phi());
0878         h_Jets_meas_no_nHCal_p->Fill(jet1_p, jet2_p);
0879         h_Jets_meas_no_nHCal_pT->Fill(jet1.pt(), jet2.pt());
0880         h_Jets_meas_no_nHCal_E->Fill(jet1.E(), jet2.E());
0881     }
0882 
0883 
0884 
0885     // measured jets vs. partons
0886 
0887     int jetid1 = -1;
0888     int jetid2 = -1;
0889 
0890     FindPartonJet(measured_jets, event[iParton_1], event[iParton_2], jetid1, jetid2);
0891 
0892     if(jetid1 >= 0)
0893     {
0894         h_Jets_meas_Partons_eta->Fill(event[iParton_1].eta(), measured_jets[jetid1].eta());
0895         h_Jets_meas_Partons_phi->Fill(event[iParton_1].phi(), measured_jets[jetid1].phi());
0896         h_Jets_meas_Partons_E->Fill(event[iParton_1].e(), measured_jets[jetid1].E());
0897 
0898         h_Jet_meas_Parton_eta1->Fill(event[iParton_1].eta(), measured_jets[jetid1].eta());
0899         h_Jet_meas_Parton_phi1->Fill(event[iParton_1].phi(), measured_jets[jetid1].phi());
0900         h_Jet_meas_Parton_E1->Fill(event[iParton_1].e(), measured_jets[jetid1].E());
0901     }
0902 
0903     if(jetid2 >= 0)
0904     {
0905         h_Jets_meas_Partons_eta->Fill(event[iParton_2].eta(), measured_jets[jetid2].eta());
0906         h_Jets_meas_Partons_phi->Fill(event[iParton_2].phi(), measured_jets[jetid2].phi());
0907         h_Jets_meas_Partons_E->Fill(event[iParton_2].e(), measured_jets[jetid2].E());
0908 
0909         h_Jet_meas_Parton_eta2->Fill(event[iParton_2].eta(), measured_jets[jetid2].eta());
0910         h_Jet_meas_Parton_phi2->Fill(event[iParton_2].phi(), measured_jets[jetid2].phi());
0911         h_Jet_meas_Parton_E2->Fill(event[iParton_2].e(), measured_jets[jetid2].E());
0912     }
0913 
0914 
0915     // events -----------------
0916 
0917 
0918     h_Event_Q2->Fill(Q2);
0919     h_Event_x->Fill(x);
0920     h_Event_y->Fill(y);
0921 
0922     // events with at least 2 jets
0923     if(inclusive_jets.size() >= 2)
0924     {
0925         h_Event_jets_Q2->Fill(Q2);
0926         h_Event_jets_x->Fill(x);
0927         h_Event_jets_y->Fill(y);
0928     }
0929 
0930     // jets
0931     if(nHCal_jets == 0)
0932     {
0933         h_Event_nHCal_0_Q2->Fill(Q2);
0934         h_Event_nHCal_0_x->Fill(x);
0935         h_Event_nHCal_0_y->Fill(y);
0936     }
0937 
0938     if(nHCal_jets == 1)
0939     {
0940         h_Event_nHCal_1_Q2->Fill(Q2);
0941         h_Event_nHCal_1_x->Fill(x);
0942         h_Event_nHCal_1_y->Fill(y);
0943     }
0944 
0945 
0946     if(nHCal_jets == 2)
0947     {
0948         h_Event_nHCal_2_Q2->Fill(Q2);
0949         h_Event_nHCal_2_x->Fill(x);
0950         h_Event_nHCal_2_y->Fill(y);
0951     }
0952 
0953     if(anyHcal_jets) // Jets in any HCal,  Both jets in any HCal
0954     {
0955         h_Event_AllHCal_Q2->Fill(Q2);
0956         h_Event_AllHCal_x->Fill(x);
0957         h_Event_AllHCal_y->Fill(y);
0958     }
0959 
0960 
0961     // measured jets
0962     if(nHCal_jets_meas == 0)
0963     {
0964         h_Event_JetMeas_nHCal_0_Q2->Fill(Q2);
0965         h_Event_JetMeas_nHCal_0_x->Fill(x);
0966         h_Event_JetMeas_nHCal_0_y->Fill(y);
0967     }
0968 
0969     if(nHCal_jets_meas == 1)
0970     {
0971         h_Event_JetMeas_nHCal_1_Q2->Fill(Q2);
0972         h_Event_JetMeas_nHCal_1_x->Fill(x);
0973         h_Event_JetMeas_nHCal_1_y->Fill(y);
0974     }
0975 
0976 
0977     if(nHCal_jets_meas == 2)
0978     {
0979         h_Event_JetMeas_nHCal_2_Q2->Fill(Q2);
0980         h_Event_JetMeas_nHCal_2_x->Fill(x);
0981         h_Event_JetMeas_nHCal_2_y->Fill(y);
0982     }
0983 
0984     if(anyHcal_jets_meas) // Jets in any HCal,  Both jets in any HCal
0985     {
0986         h_Event_JetMeas_AllHCal_Q2->Fill(Q2);
0987         h_Event_JetMeas_AllHCal_x->Fill(x);
0988         h_Event_JetMeas_AllHCal_y->Fill(y);
0989     }
0990 
0991 
0992     if(anyHcal_jets_meas_no_nHCal) // Jets in any HCal, no nHCal
0993     {
0994         h_Event_JetMeas_no_nHCal_Q2->Fill(Q2);
0995         h_Event_JetMeas_no_nHCal_x->Fill(x);
0996         h_Event_JetMeas_no_nHCal_y->Fill(y);
0997     }
0998 
0999 
1000 
1001     // tell the user what was done
1002     //  - the description of the algorithm used
1003     //  - extract the inclusive jets with pt > 5 GeV
1004     //    show the output as
1005     //      {index, rap, phi, pt}
1006     //----------------------------------------------------------
1007     //cout << "Ran " << jet_def.description() << endl;
1008 
1009     return 1;
1010 }
1011 
1012 //
1013 //  Acceptance filter
1014 //
1015 bool isInAcceptance(int i, const Event& event)
1016 {
1017     // accept all (useful for many studies)
1018     // return true;
1019     
1020     // limit to STAR TPC/BEMC/ToF acceptance + EEMC
1021     double eta = event[i].eta();
1022     if (fabs(eta) < 2)
1023         return true;
1024     else
1025         return false;
1026 }
1027 
1028 //
1029 // Cosine of angle of electron daughter in J/Psi rest frame
1030 // i.e, cos(theta)* of decay
1031 // 
1032 // Reference frame: recoil
1033 // 
1034 // Returns in range [0,1]
1035 //
1036 // Polarization: dN/dcost* = 1+alpha*cost*^2
1037 // alpha = +1 means tranverse (helicity = +-1) 
1038 //       = -1 means long. (helicity 0)
1039 //       = 0  unpolarized
1040 
1041 double costhetastar(int im, int ie, const Event& event)
1042 {
1043     double gammaFrame = event[im].e()/event[im].m(); // gamma = E/m
1044     double gammabetaFrame = event[im].pAbs()/event[im].m();  // gamma*beta = p/m
1045     
1046     double costheta = (event[im].px()*event[ie].px()+
1047                        event[im].py()*event[ie].py()+
1048                        event[im].pz()*event[ie].pz())/(event[im].pAbs()*event[ie].pAbs());
1049     double pl = event[ie].pAbs()*costheta;    // wrt J/Psi axis
1050     double pt = sin(acos(costheta))*event[ie].pAbs();   // wrt J/Psi axis
1051     double pzstar = -gammabetaFrame*event[ie].e()+gammaFrame*pl;
1052     double tanThetaStar = pt/pzstar;
1053     double thetaStar = atan(tanThetaStar);
1054     return cos(thetaStar);
1055 }
1056 
1057 int findFinalElectron(const Event& event)
1058 {
1059 
1060     int eleid = 2;
1061     bool found = false;
1062 
1063     while(found==false)
1064     {
1065         int d1 = event[eleid].daughter1();
1066         int d2 = event[eleid].daughter2();
1067 
1068 
1069         if(event[d1].id() == 11) eleid = event[eleid].daughter1();
1070         if(event[d2].id() == 11) eleid = event[eleid].daughter2();
1071 
1072         if(event[eleid].daughter1() == event[eleid].daughter2()) found = true;
1073 
1074     }
1075 
1076     return eleid;
1077 }
1078 
1079 
1080 int FillHCals(TH2F *hist, TH1F *hEne, TH1F *hEneDenom, bool &anyHcal_jets)
1081 {
1082 
1083     TH1F *hEne_Norm = (TH1F *)hEne->Clone("hEne_Norm");
1084 
1085 
1086     float integral = hEne->Integral();
1087 
1088     //cout<<"integral = "<<integral<<endl;
1089     //cout<<"Esum = "<<hEneDenom->GetBinContent(1)<<endl;
1090 
1091     hEne_Norm->Divide(hEneDenom);
1092 
1093     std::map<int, int> BinIdToCalo = { {1,-1}, {2,0}, {3,-1}, {4,1}, {5,-1}, {6,2}, {7,-1} };
1094 
1095 
1096     std::vector<pair<int, double>> input;
1097     std::vector<pair<int, double>> sorted;
1098 
1099     for (int i = 1; i <= hEne_Norm->GetNbinsX(); ++i) {
1100 
1101         pair<int, double> binPair;
1102 
1103         binPair.first = i;
1104         binPair.second = hEne_Norm->GetBinContent(i);
1105 
1106         input.push_back(binPair);
1107     }
1108 
1109     SortPairs(input, sorted);
1110 
1111     if(sorted.at(0).second > 2.0*sorted.at(1).second)
1112     {
1113         sorted.at(1) = sorted.at(0);
1114     }
1115 
1116 
1117     //cout<<"Final:"<<endl;
1118     //PrintVec(sorted);
1119 
1120     int a = BinIdToCalo.at(sorted.at(0).first);
1121     int b = BinIdToCalo.at(sorted.at(1).first);
1122 
1123     hist->Fill(a, b);
1124 
1125     if(b >= 0)
1126     {
1127         hist->Fill(a, 3);
1128     }
1129     if(a >= 0)
1130     {
1131         hist->Fill(3, b);
1132     }
1133     if(a >= 0 && b >= 0)
1134     {
1135         hist->Fill(3, 3);
1136     }
1137 
1138     a = BinIdToCalo.at(sorted.at(0).first);
1139     b = BinIdToCalo.at(sorted.at(1).first);
1140 
1141     //cout<<"a ="<<a<<endl;
1142     //cout<<"b ="<<b<<endl;
1143 
1144     if(a >= 0 && b >= 0 )
1145     {
1146         anyHcal_jets = true;
1147         //cout<<"Both jets in any HCal"<<endl;
1148     }
1149 
1150     if(a == 0 && b == 0) return 2; // 2 jets in nHCal
1151     if(a == 0 || b == 0) return 1; // 1 jet in nHCal
1152     if(a != 0 && b != 0) return 0; // 0 jet in nHCal
1153     //if(a == 3 && b == 3) return 3; // jets in any HCal
1154 
1155     return -1;
1156 
1157 }
1158 
1159 
1160 int FillHCalsJets(TH2F *hist, vector<fastjet::PseudoJet> jets)
1161 {
1162 
1163 
1164     bool anyHcal_jets = false;
1165     int nHCal_jets = 0;
1166     int a = -1;
1167     int b = -1;
1168 
1169     for (unsigned int i = 0; i < jets.size(); i++) {
1170 
1171         fastjet::PseudoJet jet = jets[i];
1172 
1173         // nHCal
1174         if(-4.14 < jet.eta() && jet.eta() < -1.18)
1175         {
1176             anyHcal_jets = true;
1177             nHCal_jets++;
1178 
1179             if(i==0) a = 0;
1180             if(i==1) b = 0;
1181         }
1182 
1183         // bHCal
1184         if(-1.1 < jet.eta() && jet.eta() < 1.1)
1185         {
1186             anyHcal_jets = true;
1187 
1188             if(i==0) a = 1;
1189             if(i==1) b = 1;
1190         }
1191 
1192         // LFHCal
1193         if(1.2 < jet.eta() && jet.eta() < 4.2)
1194         {
1195             anyHcal_jets = true;
1196 
1197             if(i==0) a = 2;
1198             if(i==1) b = 2;
1199         }
1200 
1201     }
1202 
1203     hist->Fill(a, b);
1204 
1205     if(b >= 0)
1206     {
1207         hist->Fill(a, 3);
1208     }
1209     if(a >= 0)
1210     {
1211         hist->Fill(3, b);
1212     }
1213     if(a >= 0 && b >= 0)
1214     {
1215         hist->Fill(3, 3);
1216     }
1217 
1218     //cout<<"a ="<<a<<endl;
1219     //cout<<"b ="<<b<<endl;
1220 
1221     if(a >= 0 && b >= 0 )
1222     {
1223         anyHcal_jets = true;
1224         //cout<<"Both jets in any HCal"<<endl;
1225     }
1226 
1227     if(a == 0 && b == 0) return 2; // 2 jets in nHCal
1228     if(a == 0 || b == 0) return 1; // 1 jet in nHCal
1229     if(a != 0 && b != 0) return 0; // 0 jet in nHCal
1230     //if(a == 3 && b == 3) return 3; // jets in any HCal
1231 
1232     return -1;
1233 
1234 }
1235 
1236 
1237 int FillHCalsJetsShare(TH2F *hist, vector<fastjet::PseudoJet> jets)
1238 {
1239 
1240 
1241     bool anyHcal_jets = false;
1242     int nHCal_jets = 0;
1243     int a = -1;
1244     int b = -1;
1245 
1246     bool is_nHCal = false;
1247     bool is_bHCal = false;
1248     bool is_LFHCal = false;
1249 
1250     for (unsigned int i = 0; i < jets.size(); i++) {
1251 
1252         fastjet::PseudoJet jet = jets[i];
1253 
1254         vector<PseudoJet> constituents = jet.constituents();
1255 
1256         for (int j = 0; j < constituents.size(); ++j) {
1257 
1258             fastjet::PseudoJet con = constituents[j];
1259 
1260             if(-4.14 < con.eta() && con.eta() < -1.18) is_nHCal = true;
1261             if(-1.1 < con.eta() && con.eta() < 1.1) is_bHCal = true;
1262             if(1.2 < con.eta() && con.eta() < 4.2) is_LFHCal = true;
1263         }
1264 
1265 
1266         //----------
1267 
1268         // nHCal
1269         if(is_nHCal && !is_bHCal && !is_LFHCal)
1270         {
1271             if(i==0) a = 0;
1272             if(i==1) b = 0;
1273         }
1274 
1275         // nHCal+bHCal
1276         if(is_nHCal && is_bHCal && !is_LFHCal)
1277         {
1278             if(i==0) a = 1;
1279             if(i==1) b = 1;
1280         }
1281 
1282         // bHCal
1283         if(!is_nHCal && is_bHCal && !is_LFHCal)
1284         {
1285             if(i==0) a = 2;
1286             if(i==1) b = 2;
1287         }
1288 
1289         // bHCal+LFHCal
1290         if(!is_nHCal && is_bHCal && is_LFHCal)
1291         {
1292             if(i==0) a = 3;
1293             if(i==1) b = 3;
1294         }
1295 
1296         // LFHCal
1297         if(!is_nHCal && !is_bHCal && is_LFHCal)
1298         {
1299             if(i==0) a = 4;
1300             if(i==1) b = 4;
1301         }
1302 
1303     }
1304 
1305     hist->Fill(a, b);
1306 
1307     if(b >= 0)
1308     {
1309         hist->Fill(a, 5);
1310     }
1311     if(a >= 0)
1312     {
1313         hist->Fill(5, b);
1314     }
1315     if(a >= 0 && b >= 0)
1316     {
1317         hist->Fill(5, 5);
1318     }
1319 
1320     //cout<<"a ="<<a<<endl;
1321     //cout<<"b ="<<b<<endl;
1322 
1323     if(a >= 0 && b >= 0 )
1324     {
1325         anyHcal_jets = true;
1326         //cout<<"Both jets in any HCal"<<endl;
1327     }
1328 
1329     //if(a == 0 && b == 0) return 2; // 2 jets in nHCal
1330     //if(a == 0 || b == 0) return 1; // 1 jet in nHCal
1331     //if(a != 0 && b != 0) return 0; // 0 jet in nHCal
1332     //if(a == 3 && b == 3) return 3; // jets in any HCal
1333 
1334     return -1;
1335 
1336 }
1337 
1338 
1339 void SortPairs(std::vector<pair<int, double>> &input, std::vector<pair<int, double>> &output)
1340 {
1341 
1342     //int maxId = 0;
1343     //int maxVal = 0.0;
1344 
1345     std::vector<pair<int, double>> input_copy = input;
1346 
1347     //PrintVec(input_copy);
1348 
1349     sort(input_copy.rbegin(), input_copy.rend(), compFunc);
1350 
1351     //cout<<"Sorted:"<<endl;
1352     //PrintVec(input_copy);
1353 
1354     output = input_copy;
1355 
1356 /*
1357     for (int i = 0; i < input_copy.size(); ++i) {
1358 
1359         for (int j = 0; j < input_copy.size(); ++j) {
1360 
1361             double currVal = input_copy.at(j).second;
1362 
1363             if(currVal > maxVal)
1364             {
1365                 maxId = j+1;
1366                 maxVal = currVal;
1367             }
1368         } // j
1369 
1370         pair<int, double> newPair;
1371         newPair.first = maxId;
1372         newPair.second = maxVal;
1373 
1374 
1375         output.push_back(newPair);
1376         input_copy.erase(maxId-1);
1377 
1378     } // i
1379 */
1380 }
1381 
1382 
1383 
1384 bool compFunc(pair<int, double> a, pair<int, double> b)
1385 {
1386     return (a.second < b.second);
1387 }
1388 
1389 
1390 void PrintVec(std::vector<pair<int, double>> input)
1391 {
1392 
1393     for (int i = 0; i < input.size(); ++i) {
1394 
1395         int a = input.at(i).first;
1396         float b = input.at(i).second;
1397 
1398         cout<<i<<": <"<<a<<","<<b<<">\t";
1399 
1400     }
1401     cout<<endl;
1402 
1403 }
1404 
1405 int isInHcals(double eta)
1406 {
1407 
1408     // nHCal
1409     if(-4.14 < eta && eta < -1.18)
1410     {
1411         return 1;
1412     }
1413 
1414     // bHCal
1415     if(-1.1 < eta && eta < 1.1)
1416     {
1417         return 2;
1418     }
1419 
1420     // LFHCal
1421     if(1.2 < eta && eta < 4.2)
1422     {
1423         return 3;
1424     }
1425 
1426     return 0;
1427 }
1428 
1429 
1430 bool isInTracker(double eta)
1431 {
1432     // MPGD https://indico.bnl.gov/event/20727/contributions/93067/
1433     if(-3.61 < eta && eta < 3.44)
1434     {
1435         return true;
1436     }
1437 
1438     return false;
1439 }
1440 
1441 
1442 void FindPartonJet(vector<fastjet::PseudoJet> jets, Particle parton1, Particle parton2, int &jetid1, int &jetid2)
1443 {
1444 
1445     double r1_min = 1.0;
1446     double r2_min = 1.0;
1447 
1448     int id1_min = -1;
1449     int id2_min = -1;
1450 
1451     for (int i = 0; i < jets.size(); ++i) {
1452 
1453         fastjet::PseudoJet jet = jets[i];
1454 
1455         double r1 = dR(jet, parton1);
1456         double r2 = dR(jet, parton2);
1457 
1458         if(r1<r1_min)
1459         {
1460             r1_min = r1;
1461             id1_min = i;
1462         }
1463         if(r2<r2_min)
1464         {
1465             r2_min = r2;
1466             id2_min = i;
1467         }
1468 
1469         //cout<<"r1_min = "<<r1_min<<"\tr1 = "<<r1<<endl;
1470         //cout<<"r2_min = "<<r2_min<<"\tr2 = "<<r2<<endl;
1471 
1472     } // jets
1473 
1474     jetid1 = id1_min;
1475     jetid2 = id2_min;
1476 
1477     //cout<<"jetid1 = "<<jetid1<<endl;
1478     //cout<<"jetid2 = "<<jetid2<<endl;
1479 
1480 }
1481 
1482 
1483 double dR(fastjet::PseudoJet jet, Particle parton)
1484 {
1485 
1486     double d1sq = TMath::Power(jet.eta()-parton.eta(), 2.0);
1487     d1sq += TMath::Power(jet.phi()-parton.phi(), 2.0);
1488 
1489     return TMath::Sqrt(d1sq);
1490 }
[ 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!