Pi0MassAnalysis/src/Pi0MassAnalysis.cc

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 Pi0MassAnalysis.h.
0007 //
0008 // Pi0MassAnalysis(const std::string &name = "Pi0MassAnalysis")
0009 // everything is keyed to Pi0MassAnalysis, 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 // Pi0MassAnalysis::~Pi0MassAnalysis()
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 Pi0MassAnalysis::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 Pi0MassAnalysis::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 Pi0MassAnalysis::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 Pi0MassAnalysis::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 Pi0MassAnalysis::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 Pi0MassAnalysis::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 Pi0MassAnalysis::Reset(PHCompositeNode *topNode)
0057 // not really used - it is called before the dtor is called
0058 //
0059 // void Pi0MassAnalysis::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 "Pi0MassAnalysis.h"
0065
0066 #include <calobase/RawCluster.h>
0067 #include <calobase/RawClusterContainer.h>
0068 #include <calobase/RawClusterUtility.h>
0069
0070 #include <fun4all/Fun4AllReturnCodes.h>
0071
0072 #include <globalvertex/GlobalVertex.h>
0073 #include <globalvertex/GlobalVertexMap.h>
0074
0075 #include <phool/PHCompositeNode.h>
0076 #include <phool/getClass.h>
0077
0078 #include <boost/format.hpp>
0079
0080 #include "TLorentzVector.h"
0081
0082 //____________________________________________________________________________..
0083 Pi0MassAnalysis::Pi0MassAnalysis(const std::string &name)
0084   : SubsysReco(name)
0085 {
0086   std::cout << "Pi0MassAnalysis::Pi0MassAnalysis(const std::string &name) Calling ctor" << std::endl;
0087   _foutname = name;
0088 }
0089
0090 //____________________________________________________________________________..
0091 Pi0MassAnalysis::~Pi0MassAnalysis()
0092 {
0093   std::cout << "Pi0MassAnalysis::~Pi0MassAnalysis() Calling dtor" << std::endl;
0094 }
0095
0096 //____________________________________________________________________________..
0097 int Pi0MassAnalysis::Init(PHCompositeNode *topNode)
0098 {
0099   std::cout << "Pi0MassAnalysis::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0100
0101   counter = 0;
0102
0103   _f = new TFile(_foutname.c_str(), "RECREATE");
0104
0105   _tree = new TTree("ttree", "");
0106
0107   _tree->Branch("nclus", &nclus, "nclus/I");
0108   _tree->Branch("cluster_energy", cluster_energy, "cluster_energy[nclus]/F");
0109   _tree->Branch("cluster_eta", cluster_eta, "cluster_eta[nclus]/F");
0110   _tree->Branch("cluster_phi", cluster_phi, "cluster_phi[nclus]/F");
0111   _tree->Branch("cluster_prob", cluster_prob, "cluster_prob[nclus]/F");
0112   _tree->Branch("cluster_chi2", cluster_chi2, "cluster_chi2[nclus]/F");
0113
0114   _tree->Branch("npi0", &npi0, "npi0/I");
0115   _tree->Branch("pi0cand_pt", pi0cand_pt, "pi0cand_pt[npi0]/F");
0116   _tree->Branch("pi0cand_eta", pi0cand_eta, "pi0cand_eta[npi0]/F");
0117   _tree->Branch("pi0cand_phi", pi0cand_phi, "pi0cand_phi[npi0]/F");
0118   _tree->Branch("pi0cand_mass", pi0cand_mass, "pi0cand_mass[npi0]/F");
0119
0120   pi0MassHist = new TH1F("pi0mass", ";m_{#gamma#gamma} [GeV];Entries", 50, 0.0, 0.3);
0121
0122   for (int i = 0; i < 24; i++)
0123   {
0124     boost::format formatting("pi0mass %.1lf < eta < %.1lf");
0125     std::string histName = (formatting % (-1.2 + i * 0.1) % (-1.1 + i * 0.1)).str();
0126     pi0MassHistEtaDep[i] = new TH1F(histName.c_str(), ";m_{#gamma#gamma} [GeV];Entries", 50, 0.0, 0.3);
0127   }
0128
0129   return Fun4AllReturnCodes::EVENT_OK;
0130 }
0131
0132 //____________________________________________________________________________..
0133 int Pi0MassAnalysis::InitRun(PHCompositeNode *topNode)
0134 {
0135   std::cout << "Pi0MassAnalysis::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0136   return Fun4AllReturnCodes::EVENT_OK;
0137 }
0138
0139 //____________________________________________________________________________..
0140 int Pi0MassAnalysis::process_event(PHCompositeNode *topNode)
0141 {
0142   nclus = 0;
0143
0144   // Get Vertex
0145   float vx = 0;
0146   float vy = 0;
0147   float vz = 0;
0148   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0149   if (vertexmap)
0150   {
0151     if (!vertexmap->empty())
0152     {
0153       GlobalVertex *vtx = (vertexmap->begin()->second);
0154       vx = vtx->get_x();
0155       vy = vtx->get_y();
0156       vz = vtx->get_z();
0157     }
0158   }
0159
0160   CLHEP::Hep3Vector vertex(vx, vy, vz);
0161
0162   RawClusterContainer *clusterEM = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
0163   if (clusterEM)
0164   {
0165     RawClusterContainer::Range begin_end = clusterEM->getClusters();
0166     for (RawClusterContainer::Iterator iter = begin_end.first; iter != begin_end.second; ++iter)
0167     {
0168       RawCluster *clust = iter->second;
0169
0170       if (clust->get_energy() > minClusterEnergy && clust->get_prob() > photonClusterProbability)
0171       {
0172         // CLHEP::Hep3Vector origin(0, 0, 0);
0173         CLHEP::Hep3Vector cluster_vector = RawClusterUtility::GetECoreVec(*clust, vertex);
0174
0175         cluster_energy[nclus] = clust->get_energy();
0176         cluster_eta[nclus] = cluster_vector.pseudoRapidity();
0177         cluster_phi[nclus] = clust->get_phi();
0178         cluster_prob[nclus] = clust->get_prob();
0179         cluster_chi2[nclus] = clust->get_chi2();
0180
0181         nclus++;
0182       }
0183     }
0184   }
0185
0186   npi0 = 0;
0187
0188   for (int i = 0; i < nclus; i++)
0189   {
0190     TLorentzVector v1;
0191     v1.SetPtEtaPhiM(cluster_energy[i] / cosh(cluster_eta[i]), cluster_eta[i], cluster_phi[i], 0.0);
0192     for (int j = i + 1; j < nclus; j++)
0193     {
0194       float alpha = fabs(cluster_energy[i] - cluster_energy[j]) / (cluster_energy[i] + cluster_energy[j]);
0195
0196       if (alpha > 0.5)
0197         continue;
0198
0199       TLorentzVector v2;
0200       v2.SetPtEtaPhiM(cluster_energy[j] / cosh(cluster_eta[j]), cluster_eta[j], cluster_phi[j], 0.0);
0201
0202       if (v1.DeltaR(v2) > 0.8)
0203         continue;
0204
0205       TLorentzVector res;
0206       res = v1 + v2;
0207
0208       pi0cand_pt[npi0] = res.Pt();
0209       pi0cand_eta[npi0] = res.Eta();
0210       pi0cand_phi[npi0] = res.Phi();
0211       pi0cand_mass[npi0] = res.M();
0212       npi0++;
0213
0214       pi0MassHist->Fill(res.M());
0215
0216       int eta_index = floor((res.Eta() + 1.2) * 10);
0217
0218       if (eta_index >= 0 && eta_index < 24)
0219       {
0220         pi0MassHistEtaDep[eta_index]->Fill(res.M());
0221       }
0222     }
0223   }
0224
0225   _tree->Fill();
0226
0227   counter++;
0228
0229   if (counter % 100 == 0)
0230     std::cout << counter << " events are processed" << std::endl;
0231
0232   return Fun4AllReturnCodes::EVENT_OK;
0233 }
0234
0235 //____________________________________________________________________________..
0236 int Pi0MassAnalysis::ResetEvent(PHCompositeNode *topNode)
0237 {
0238   return Fun4AllReturnCodes::EVENT_OK;
0239 }
0240
0241 //____________________________________________________________________________..
0242 int Pi0MassAnalysis::EndRun(const int runnumber)
0243 {
0244   std::cout << "Pi0MassAnalysis::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0245   return Fun4AllReturnCodes::EVENT_OK;
0246 }
0247
0248 //____________________________________________________________________________..
0249 int Pi0MassAnalysis::End(PHCompositeNode *topNode)
0250 {
0251   std::cout << "Pi0MassAnalysis::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0252
0253   _f->Write();
0254   _f->Close();
0255
0256   return Fun4AllReturnCodes::EVENT_OK;
0257 }
0258
0259 //____________________________________________________________________________..
0260 int Pi0MassAnalysis::Reset(PHCompositeNode *topNode)
0261 {
0262   std::cout << "Pi0MassAnalysis::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0263   return Fun4AllReturnCodes::EVENT_OK;
0264 }
0265
0266 //____________________________________________________________________________..
0267 void Pi0MassAnalysis::Print(const std::string &what) const
0268 {
0269   std::cout << "Pi0MassAnalysis::Print(const std::string &what) const Printing info for " << what << std::endl;
0270 }