calotriggeremulator/src/SmearTruthJets.cc

0001 #include "SmearTruthJets.h"
0002 //for emc clusters
0003 #include <TH1.h>
0004 #include <TH2.h>
0005 #include <TFile.h>
0006 #include <TVector3.h>
0007 #include <TH2D.h>
0008 #include <TF1.h>
0009 #include <TH1D.h>
0010 #include <TMath.h>
0011 #include <jetbase/Jet.h>
0012 #include <jetbase/JetContainer.h>
0013 #include <jetbase/JetContainerv1.h>
0014 #include <jetbase/Jetv2.h>
0015 #include <vector>
0016
0017 #include <fun4all/Fun4AllReturnCodes.h>
0018 #include <phool/PHCompositeNode.h>
0019 #include <phool/PHIODataNode.h>
0020 #include <phool/PHNode.h>
0021 #include <phool/PHNodeIterator.h>
0022 #include <phool/PHObject.h>
0023 #include <phool/getClass.h>
0024 // G4Cells includes
0025
0026 #include <iostream>
0027
0028 #include <map>
0029
0030 //____________________________________________________________________________..
0031 SmearTruthJets::SmearTruthJets(const std::string &name):
0032   SubsysReco(name)
0033
0034 {
0035 }
0036
0037 //____________________________________________________________________________..
0038 SmearTruthJets::~SmearTruthJets()
0039 {
0040
0041 }
0042
0043 //____________________________________________________________________________..
0044 int SmearTruthJets::Init(PHCompositeNode *topNode)
0045 {
0046   TFile *f = new TFile("/sphenix/user/dlis/Projects/CaloTriggerEmulator/analysis/calotriggeremulator/src/jesjer.root","r");
0047   h_jes = (TH1D*) f->Get("h_jes")->Clone();
0048   h_jer = (TH1D*) f->Get("h_jer")->Clone();
0049   if (!h_jer)
0050     {
0051       std::cout << "NOTHING" << std::endl;
0052     }
0053   f->Close();
0054
0055   smearfunction = new TF1("smearfunction","gaus", 0, 2);
0056   return Fun4AllReturnCodes::EVENT_OK;
0057 }
0058
0059 //____________________________________________________________________________..
0060 int SmearTruthJets::InitRun(PHCompositeNode *topNode)
0061 {
0062
0063   PHNodeIterator iter(topNode);
0064
0065   // Looking for the DST node
0066   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0067   if (!dstNode)
0068     {
0069       std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0070       return Fun4AllReturnCodes::ABORTRUN;
0071     }
0072
0073   // Create the AntiKt node if required
0074   PHCompositeNode *AlgoNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "ANTIKT"));
0075   if (!AlgoNode)
0076     {
0077       AlgoNode = new PHCompositeNode("ANTIKT");
0078       dstNode->addNode(AlgoNode);
0079     }
0080
0081   // Create the Input node if required
0082   PHCompositeNode *InputNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "TRUTH"));
0083   if (!InputNode)
0084     {
0085       InputNode = new PHCompositeNode("TRUTH");
0086       AlgoNode->addNode(InputNode);
0087     }
0088
0089   JetContainer *jetconn = findNode::getClass<JetContainer>(topNode, "AntiKt_TruthSmear_r04");
0090   if (!jetconn)
0091     {
0092       jetconn = new JetContainerv1();
0093       PHIODataNode<PHObject> *JetContainerNode = new PHIODataNode<PHObject>(jetconn, "AntiKt_TruthSmear_r04", "PHObject");
0094       InputNode->addNode(JetContainerNode);
0095     }
0096   return Fun4AllReturnCodes::EVENT_OK;
0097 }
0098
0099 int SmearTruthJets::process_event(PHCompositeNode *topNode)
0100 {
0101   std::cout << " here"<<std::endl;
0102
0103   std::string recoJetName = "AntiKt_Truth_r04";
0104   JetContainer *jets_truth_4 = findNode::getClass<JetContainer>(topNode, recoJetName);
0105   std::cout << " here"<<std::endl;
0106   recoJetName = "AntiKt_TruthSmear_r04";
0107   JetContainer *jets_truth_smear_4 = findNode::getClass<JetContainer>(topNode, recoJetName);
0108   std::cout << " here"<<std::endl;
0109   jets_truth_smear_4 = jets_truth_4;
0110   std::cout << " here"<<std::endl;
0111   for (auto jet : *jets_truth_smear_4)
0112     {
0113       std::cout << "JET"<<std::endl;
0114       std::cout << "Start: " << jet->get_pt() << " " << jet->get_eta() << " " << jet->get_phi() << std::endl;
0115       float pt = jet->get_pt();
0116       int ptbin = floor(pt/5) - 2;
0117       if (ptbin < 1) ptbin = 1;
0118       if (ptbin > 9) ptbin = 9;
0119       float jes = h_jes->GetBinContent(ptbin);
0120       float jer = h_jer->GetBinContent(ptbin);
0121
0122       smearfunction->SetParameters(1, jes, jer);
0123
0124       double x = smearfunction->GetRandom();
0125       std::cout << jes << " " << jer << " -- " << x<<std::endl;
0126       pt *= x;
0127       TVector3 v;
0128       v.SetPtEtaPhi(pt, jet->get_eta(), jet->get_phi());
0129       jet->set_px(v.X());
0130       jet->set_py(v.Y());
0131       jet->set_pz(v.Z());
0132       std::cout << "End : " << jet->get_pt() << " " << jet->get_eta() << " " << jet->get_phi() << std::endl;
0133     }
0134
0135   return Fun4AllReturnCodes::EVENT_OK;
0136 }
0137
0138
0139 //____________________________________________________________________________..
0140 int SmearTruthJets::End(PHCompositeNode *topNode)
0141 {
0142   if (Verbosity() > 0)
0143     {
0144       std::cout << "SmearTruthJets::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0145     }
0146   delete smearfunction;
0147   delete h_jer;
0148   delete h_jes;
0149   return Fun4AllReturnCodes::EVENT_OK;
0150 }
0151
0152 //____________________________________________________________________________..
0153 int SmearTruthJets::Reset(PHCompositeNode *topNode)
0154 {
0155   if (Verbosity() > 0)
0156     {
0157       std::cout << "SmearTruthJets::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0158     }
0159   return Fun4AllReturnCodes::EVENT_OK;
0160 }