mainAnalysis/src/AnNeutralMeson.cc

0001 #include "AnNeutralMeson.h"
0002
0003 // Spin DB
0004 #include <uspin/SpinDBContent.h>
0005 #include <uspin/SpinDBOutput.h>
0006
0007 // Tower includes
0008 #include <calobase/RawCluster.h>
0009 #include <calobase/RawClusterContainer.h>
0010 #include <calobase/RawClusterUtility.h>
0011 #include <calobase/RawTower.h>
0012 #include <calobase/RawTowerContainer.h>
0013 #include <calobase/RawTowerGeom.h>
0014 #include <calobase/RawTowerGeomContainer.h>
0015 #include <calobase/TowerInfo.h>
0016 #include <calobase/TowerInfoContainer.h>
0017 #include <calobase/TowerInfoDefs.h>
0018
0019 #include <fun4all/Fun4AllHistoManager.h>
0020 #include <fun4all/Fun4AllReturnCodes.h>
0021
0022 #include <phool/getClass.h>
0023 #include <phool/PHCompositeNode.h>
0024
0025 #include <globalvertex/GlobalVertex.h>
0026 #include <globalvertex/MbdVertex.h>
0027
0028 // MBD
0029 #include <mbd/MbdGeom.h>
0030 #include <mbd/MbdPmtContainerV1.h>
0031 #include <mbd/MbdPmtHit.h>
0032
0033 #include <Math/Vector4D.h>
0034 #include <TFile.h>
0035 #include <TH1.h>
0036 #include <TH2.h>
0037
0038 #include <Event/Event.h>
0039 #include <Event/packet.h>
0040 #include <cassert>
0041 #include <fstream>
0042 #include <sstream>
0043 #include <string>
0044
0045 // To check virtual and resident memory usage
0046 //#include "monitorMemory.h"
0047
0048 #include <algorithm>
0049 #include <numeric>
0050 #include <fstream>
0051 #include <iostream>
0052 #include <utility>
0053 #include "TCanvas.h"
0054 #include "TF1.h"
0055 #include "TFile.h"
0056 #include "TH1F.h"
0057 #include "TH2F.h"
0058 #include "TH3F.h"
0059 #include "TRandom3.h"
0060
0061 #include <cdbobjects/CDBHistos.h>  // for CDBHistos
0062 #include <cdbobjects/CDBTTree.h>   // for CDBTTree
0063
0064 #include <ffamodules/CDBInterface.h>
0065 #include <phool/recoConsts.h>
0066
0067 #include <g4main/PHG4TruthInfoContainer.h>
0068
0069 using namespace std;
0070
0071 AnNeutralMeson::AnNeutralMeson(const std::string &name, const int &runnumber,
0072                                const std::string &filename, const std::string& foldername)
0073   : SubsysReco(name)
0074   , _runnumber(runnumber)
0075   , outfilename(filename)
0076   , outfolder(foldername)
0077 {
0078   _eventcounter = 0;
0079 }
0080
0081 AnNeutralMeson::~AnNeutralMeson()
0082 {
0083   // delete hm;
0084   delete outfile_histograms;
0085 }
0086
0087 int AnNeutralMeson::Init(PHCompositeNode *)
0088 {
0089   // create and register your histos (all types) here
0090   outfile_histograms =
0091       new TFile((outfolder + "/OUTHIST_" + outfilename).c_str(),
0092                 "RECREATE");
0093
0094   // Event QA
0095   h_event_vtx_z =
0096       new TH1F(
0097           "h_event_vtx_z",
0098           ";counts; vtx_z [cm]", 200, -200, 200);
0099
0100   // cluster QA
0101   h_clusE = new TH1F(
0102       "h_clusE",
0103       "; E [GeV]; counts", 300, 0, 15);
0104   h_clus_eta = new TH1F(
0105       "h_clus_eta",
0106       ";eta;counts", 200, -1.2, 1.2);
0107   h_clus_phi =
0108       new TH1F(
0109           "h_clus_phi",
0110           ";#phi; counts", 200, -TMath::Pi(), TMath::Pi());
0111   h_clus_eta_phi = new TH2F(
0112       "h_clus_eta_phi",
0113       ";#eta;#phi;counts", 200, -2, 2,
0114       160, -TMath::Pi(), TMath::Pi());
0115   h_clus_eta_E =
0116       new TH2F(
0117           "h_clus_eta_E",
0118           ";#eta;E;counts", 200, -2, 2, 150, 0, 15);
0119   h_clus_eta_vtxz = new TH2F(
0120       "h_clus_eta_vtxz",
0121       ";#eta;vertex z [cm];counts",
0122       200, -2, 2, 200, -200, 200);
0123   h_clus_pt = new TH1F(
0124       "h_clus_pt",
0125       ";p_{T};counts", 100, 0, 10);
0126   h_clus_chisq = new TH1F(
0127       "h_clus_chisq",
0128       ";#chi^{2};counts", 100, 0, 10);
0129
0130   return 0;
0131 }
0132
0133 int AnNeutralMeson::InitRun(PHCompositeNode *topNode)
0134 {
0135   syncobject = findNode::getClass<SyncObject>(topNode, syncdefs::SYNCNODENAME);
0136   if (!syncobject)
0137   {
0138     // Not used here but sync object should be present to keep the event number
0139     // in the output µDST
0140     std::cout << "syncobject not found" << std::endl;
0141   }
0142
0143   gl1packet = findNode::getClass<Gl1Packet>(topNode, "14001");
0144   if (!gl1packet)
0145   {
0146     std::cerr << "AnNeutralMeson Gl1Packet node is missing" << std::endl;
0147     return Fun4AllReturnCodes::ABORTRUN;
0148   }
0149
0150   vertexmap =
0151     findNode::getClass<GlobalVertexMap>(topNode,
0152                                         "GlobalVertexMap");
0153   if (!vertexmap)
0154   {
0155       std::cout << "AnNeutralMeson GlobalVertexMap node is missing"
0156                 << std::endl;
0157       return Fun4AllReturnCodes::ABORTRUN;
0158   }
0159
0160   clusterContainer =
0161     findNode::getClass<RawClusterContainer>(topNode,
0162                                             "CLUSTERINFO_CEMC");
0163   if (!clusterContainer)
0164     {
0165       std::cout << PHWHERE << "AnNeutralMeson - Fatal Error - "
0166         "CLUSTER_CEMC node is missing. "
0167                 << std::endl;
0168       return Fun4AllReturnCodes::ABORTRUN;
0169   }
0170
0171   try
0172   {
0173     CreateNodes(topNode);
0174   }
0175   catch (std::exception &e)
0176   {
0177     std::cout << PHWHERE << ": " << e.what() << std::endl;
0178     throw;
0179   }
0180
0181   return 0;
0182 }
0183
0184 int AnNeutralMeson::process_event(PHCompositeNode *topNode)
0185 {
0186   _eventcounter++;
0187
0188   return process_towers(topNode);
0189
0190   return 0;
0191 }
0192
0193 int AnNeutralMeson::process_towers(PHCompositeNode */*topNode*/)
0194 {
0195   // Show progression
0196   if ((_eventcounter % 1000) == 0)
0197     std::cout << _eventcounter << std::endl;
0198
0199   //-----------------------check MBD
0200   // trigger----------------------------------------//
0201
0202   // see Joe Mead's GL1/GTM manual (Firmware version 51)
0203   // Live counter increments only when the "busy" signal is inactive
0204   live_trigger = gl1packet->getLiveVector();
0205   // Scaledown counter only increments when the number of triggers reaches the
0206   // scaledown vector (scaledown is the ratio between live and scaled values)
0207   scaled_trigger = gl1packet->getScaledVector();
0208
0209   // 0x1U is 1 in hexadecimal notation. U means it is unsigned
0210   if ((live_trigger >> 10 & 0x1U) != 0x1U)
0211   {
0212     // Ensure that the 10th bit is equal to 1 (MBD N&S >= 1)
0213     return Fun4AllReturnCodes::ABORTEVENT;
0214   }
0215   _eventcounter_selection1++;
0216
0217   //-----------------------spin
0218   // information----------------------------------------//
0219
0220   // Only keep the bunchnumber at this point.
0221   // The spin direction and average polarization can be
0222   // deduced later at the run level (or at the fill level).
0223   bunchnumber = gl1packet->getBunchNumber();
0224
0225   //-----------------------get vertex----------------------------------------//
0226
0227   vtx_z = 0;
0228
0229   if (vertexmap && !vertexmap->empty())
0230   {
0231     GlobalVertex *vtx = vertexmap->begin()->second;
0232     if (vtx)
0233     {
0234       vtx_z = vtx->get_z();
0235       h_event_vtx_z->Fill(vtx_z);
0236     }
0237     else
0238     {
0239       return Fun4AllReturnCodes::ABORTEVENT;
0240     }
0241   }
0242   else
0243   {
0244     return Fun4AllReturnCodes::ABORTEVENT;
0245   }
0246
0247   _eventcounter_selection2++;
0248
0249   // Make sure ClusterSmallInfo container is empty before filling it with new clusters
0250   _smallclusters->Reset();
0251   _smallclusters->set_live_trigger(live_trigger);
0252   _smallclusters->set_scaled_trigger(scaled_trigger);
0253   _smallclusters->set_bunchnumber(bunchnumber);
0254
0255   // Select good clusters
0256   int num_good_clusters = 0;
0257   RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
0258   RawClusterContainer::ConstIterator clusterIter;
0259   for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second;
0260        clusterIter++)
0261   {
0262     RawCluster *recoCluster = clusterIter->second;
0263
0264     CLHEP::Hep3Vector vertex(0, 0, vtx_z);
0265     CLHEP::Hep3Vector E_vec_cluster =
0266         RawClusterUtility::GetECoreVec(*recoCluster, vertex);
0267     CLHEP::Hep3Vector pos_vec_cluster = recoCluster->get_position() - vertex;
0268
0269     float clusE = E_vec_cluster.mag();
0270     float clus_eta = E_vec_cluster.pseudoRapidity();
0271     float clus_phi = E_vec_cluster.phi();
0272     float clus_pt = E_vec_cluster.perp();
0273     float clus_chisq = recoCluster->get_chi2();
0274
0275     h_clusE->Fill(clusE);
0276     h_clus_eta->Fill(clus_eta);
0277     h_clus_phi->Fill(clus_phi);
0278     h_clus_eta_phi->Fill(clus_eta, clus_phi);
0279     h_clus_eta_E->Fill(clus_eta, clusE);
0280     h_clus_eta_vtxz->Fill(clus_eta, vtx_z);
0281     h_clus_pt->Fill(clus_pt);
0282     h_clus_chisq->Fill(clus_chisq);
0283
0284     _clustercounter++; // Cluster counter before cuts
0285
0286     if (clusE < clus_E_cut)
0287       continue;
0288     if (clus_chisq > clus_chisq_cut)
0289       continue;
0290
0291     _clustercounter_selection++; // Cluster counter after cuts
0292     num_good_clusters++;
0293
0294     // Don't store events with more than 50 clusters
0295     // Too noisy
0296     if (!(_smallclusters->add_cluster(clus_eta,
0297                                       clus_phi,
0298                                       clusE,
0299                                       recoCluster->get_energy(),
0300                                       clus_chisq)))
0301     {
0302       return Fun4AllReturnCodes::ABORTEVENT;
0303     }
0304   }  // cluster loop
0305
0306   // Compress the ClusterSmallInfoContainer object
0307   _smallclusters->compress();
0308
0309   // Don't store events with no clusters
0310   if (num_good_clusters == 0)
0311   {
0312     return Fun4AllReturnCodes::ABORTEVENT;
0313   }
0314
0315   // Store any other events
0316   return Fun4AllReturnCodes::EVENT_OK;
0317 }
0318
0319 int AnNeutralMeson::End(PHCompositeNode * /*topNode*/)
0320 {
0321   // Save QA histograms (before good cluster selection)
0322   outfile_histograms->cd();
0323   outfile_histograms->Write();
0324   outfile_histograms->Close();
0325   delete outfile_histograms;
0326
0327   // Count events and clusters
0328   std::cout
0329       << "selection (all > MBDTriggerLive > GlobalVertex): "
0330       << _eventcounter << " > " << _eventcounter_selection1 << " > "
0331       << _eventcounter_selection2 << std::endl;
0332
0333   std::cout << "selection clusters (all > (E>0.5, chi2<1000)): "
0334             << _clustercounter << " > " << _clustercounter_selection
0335             << std::endl;
0336
0337   return 0;
0338 }
0339
0340 void AnNeutralMeson::CreateNodes(PHCompositeNode *topNode)
0341 {
0342   PHNodeIterator iter(topNode);
0343
0344   // Grab the CEMC node
0345   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0346   if (!dstNode)
0347   {
0348     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0349     throw std::runtime_error("Failed to find DST node in EmcRawTowerBuilder::CreateNodes");
0350   }
0351
0352   // Get the _det_name subnode
0353   PHCompositeNode *cemcNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "CEMC"));
0354
0355   // Check that it is there
0356   if (!cemcNode)
0357   {
0358     cemcNode = new PHCompositeNode("CEMC");
0359     dstNode->addNode(cemcNode);
0360   }
0361   ClusterSmallInfoNodeName = "CLUSTER_SMALLINFO_CEMC";
0362
0363   _smallclusters = findNode::getClass<ClusterSmallInfoContainer>(dstNode, ClusterSmallInfoNodeName);
0364
0365   if (!_smallclusters)
0366   {
0367     _smallclusters = new ClusterSmallInfoContainer();
0368   }
0369
0370   PHIODataNode<PHObject> *clusterNode = new PHIODataNode<PHObject>(_smallclusters, ClusterSmallInfoNodeName, "PHObject");
0371
0372   cemcNode->addNode(clusterNode);
0373 }