neutralMesonTSSA/src/neutralMesonTSSA.h

0001 // Tell emacs that this is a C++ source
0002 //  -*- C++ -*-.
0003 #ifndef NEUTRALMESONTSSA_H
0004 #define NEUTRALMESONTSSA_H
0005
0006 #include <fun4all/SubsysReco.h>
0007
0008 #include <string>
0009 #include <vector>
0010 #include <utility>
0011
0012 class PHCompositeNode;
0013 class RunHeader;
0014 class Gl1Packet;
0015 class RawClusterContainer;
0016 class PHG4TruthInfoContainer;
0017 class GlobalVertex;
0018 class PHG4VtxPoint;
0019 class TFile;
0020 class TTree;
0021 class TH1;
0022 class TH2;
0023 class BinnedHistSet;
0024
0025 class Diphoton
0026 {
0027   public:
0028     float mass, E, eta, phi, pT, xF, vtxz;
0029     int clus1index, clus2index;
0030     /* Diphoton() {} */
0031     /* ~Diphoton() {} */
0032 };
0033
0034 class PhiHists
0035 {
0036   public:
0037     BinnedHistSet *phi_pT, *phi_pT_blue_up, *phi_pT_blue_down, *phi_pT_yellow_up, *phi_pT_yellow_down;
0038     BinnedHistSet *phi_xF, *phi_xF_blue_up, *phi_xF_blue_down, *phi_xF_yellow_up, *phi_xF_yellow_down;
0039     BinnedHistSet *phi_eta, *phi_eta_blue_up, *phi_eta_blue_down, *phi_eta_yellow_up, *phi_eta_yellow_down;
0040     BinnedHistSet *phi_vtxz, *phi_vtxz_blue_up, *phi_vtxz_blue_down, *phi_vtxz_yellow_up, *phi_vtxz_yellow_down;
0041     /* PhiHists() {} */
0042     /* ~PhiHists() {} */
0043 };
0044
0045 class neutralMesonTSSA : public SubsysReco
0046 {
0047  public:
0048
0049   neutralMesonTSSA(const std::string &name = "neutralMesonTSSA", std::string histname = "neutralMesonTSSA_hists.root", std::string treename = "neutralMesonTSSA_trees.root", bool isMC = false);
0050
0051   ~neutralMesonTSSA() override;
0052
0053   /** Called during initialization.
0054       Typically this is where you can book histograms, and e.g.
0055       register them to Fun4AllServer (so they can be output to file
0056       using Fun4AllServer::dumpHistos() method).
0057    */
0058   int Init(PHCompositeNode *topNode) override;
0059
0060   /** Called for first event when run number is known.
0061       Typically this is where you may want to fetch data from
0062       database, because you know the run number. A place
0063       to book histograms which have to know the run number.
0064    */
0065   int InitRun(PHCompositeNode *topNode) override;
0066
0067   /** Called for each event.
0068       This is where you do the real work.
0069    */
0070   int process_event(PHCompositeNode *topNode) override;
0071
0072   /// Clean up internals after each event.
0073   int ResetEvent(PHCompositeNode *topNode) override;
0074
0075   /// Called at the end of each run.
0076   int EndRun(const int runnumber) override;
0077
0078   /// Called at the end of all processing.
0079   int End(PHCompositeNode *topNode) override;
0080
0081   /// Reset
0082   int Reset(PHCompositeNode * /*topNode*/) override;
0083
0084   void Print(const std::string &what = "ALL") const override;
0085
0086   float get_min_clusterE();
0087   void set_min_clusterE(float Emin);
0088   float get_max_clusterChi2();
0089   void set_max_clusterChi2(float Chi2max);
0090   bool InRange(float mass, std::pair<float, float> range);
0091
0092   // In Init()
0093   void MakePhiHists(std::string which); // which = pi0, eta, pi0bkgr, etabkgr
0094   void MakeAllHists();
0095   void MakeVectors();
0096   // In InitRun()
0097   void GetRunNum();
0098   int GetSpinInfo(); // spin patterns, beam polarization and crossing shift for this run
0099   /* void CountLumi(); */
0100   // In process_event()
0101   void GetTrigger();
0102   void GetBunchNum();
0103   void GetSpins(); // blue and yellow spins for this event
0104   void FindGoodClusters();
0105   /* void FillClusterTree(); */
0106   void FindDiphotons();
0107   /* void FillDiphotonTree(); */
0108   /* void FillPhiHists(std::string which, int index); // which = pi0, eta, pi0bkgr, etabkgr */
0109   /* void FillAllPhiHists(); */
0110   // In ResetEvent()
0111   void ClearVectors();
0112   // In End()
0113   void DeleteStuff();
0114   // Other info
0115   void PrintTrigger();
0116
0117  private:
0118   bool isMonteCarlo;
0119   std::string outfilename_hists;
0120   std::string outfilename_trees;
0121   TFile* outfile_hists = nullptr;
0122   TFile* outfile_trees = nullptr;
0123   TTree* tree_event = nullptr;
0124   TTree* tree_clusters = nullptr;
0125   TTree* tree_diphotons = nullptr;
0126
0127   // data containers
0128   RawClusterContainer* m_clusterContainer = nullptr;
0129   PHG4TruthInfoContainer* m_truthInfo = nullptr;
0130   GlobalVertex* gVtx = nullptr;
0131   PHG4VtxPoint* mcVtx = nullptr;
0132   RunHeader* runHeader = nullptr;
0133   Gl1Packet* gl1Packet = nullptr;
0134
0135   // spin info
0136   static const int NBUNCHES = 120;
0137   int runNum = 0;
0138   int spinPatternBlue[NBUNCHES] = {0};
0139   int spinPatternYellow[NBUNCHES] = {0};
0140   int bunchNum = 0;
0141   int crossingShift = 0;
0142   int sphenixBunch = 0;
0143   int bspin = 0;
0144   int yspin = 0;
0145   long long gl1pScalers[NBUNCHES] = {0};
0146   float lumiUpBlue = 0;
0147   float lumiDownBlue = 0;
0148   /* float relLumiBlue; */
0149   float polBlue = 0;
0150   float lumiUpYellow = 0;
0151   float lumiDownYellow = 0;
0152   /* float relLumiYellow; */
0153   float polYellow = 0;
0154   float crossingAngle = -999.9;
0155   float crossingAngleIntended = -999.9;
0156
0157   // trigger & event-level info
0158   bool mbdNtrigger = false;
0159   bool mbdStrigger = false;
0160   bool mbdtrigger_live = false;
0161   bool photontrigger_live = false;
0162   bool mbdtrigger_scaled = false;
0163   bool photontrigger_scaled = false;
0164   bool mbdvertex = false;
0165   bool globalvertex = false;
0166   TH1* h_nEvents = nullptr;
0167   float vtxz = -9999999.9;
0168   TH1* h_vtxz = nullptr;
0169   TH1* h_crossingAngle = nullptr;
0170
0171   // tower-level info
0172   TH2* h_towerEta_Phi_500MeV = nullptr;
0173   TH2* h_towerEta_Phi_800MeV = nullptr;
0174
0175   // cluster distributions
0176   TH1* h_nAllClusters = nullptr;
0177   TH1* h_nGoodClusters = nullptr;
0178   TH1* h_clusterE = nullptr;
0179   TH1* h_clusterEta = nullptr;
0180   TH2* h_clusterVtxz_Eta = nullptr;
0181   TH1* h_clusterPhi = nullptr;
0182   TH2* h_clusterEta_Phi = nullptr;
0183   TH1* h_clusterpT = nullptr;
0184   TH1* h_clusterxF = nullptr;
0185   TH2* h_clusterpT_xF = nullptr;
0186   TH1* h_clusterChi2 = nullptr;
0187   TH1* h_clusterChi2zoomed = nullptr;
0188   TH1* h_mesonClusterChi2 = nullptr;
0189   TH1* h_goodClusterE = nullptr;
0190   TH1* h_goodClusterEta = nullptr;
0191   TH2* h_goodClusterVtxz_Eta = nullptr;
0192   TH1* h_goodClusterPhi = nullptr;
0193   TH2* h_goodClusterEta_Phi = nullptr;
0194   TH1* h_goodClusterpT = nullptr;
0195   TH1* h_goodClusterxF = nullptr;
0196   TH2* h_goodClusterpT_xF = nullptr;
0197
0198   // diphoton distributions
0199   TH1* h_nDiphotons = nullptr;
0200   TH1* h_nRecoPi0s = nullptr;
0201   TH1* h_nRecoEtas = nullptr;
0202   TH1* h_diphotonE = nullptr;
0203   TH1* h_diphotonMass = nullptr;
0204   TH1* h_diphotonEta = nullptr;
0205   TH2* h_diphotonVtxz_Eta = nullptr;
0206   TH1* h_diphotonPhi = nullptr;
0207   TH2* h_diphotonEta_Phi = nullptr;
0208   TH1* h_diphotonpT = nullptr;
0209   TH1* h_diphotonxF = nullptr;
0210   TH2* h_diphotonpT_xF = nullptr;
0211   TH1* h_diphotonAsym = nullptr;
0212   TH1* h_diphotonDeltaR = nullptr;
0213   TH2* h_diphotonAsym_DeltaR = nullptr;
0214   BinnedHistSet* bhs_diphotonMass_pT = nullptr;
0215   BinnedHistSet* bhs_diphotonMass_xF = nullptr;
0216
0217   // clusters and cuts
0218   float min_clusterE = 0.5;
0219   float max_clusterE = 30.0;
0220   float max_clusterChi2 = 1000.0;
0221   int nAllClusters = 0;
0222   int nGoodClusters = 0;
0223   std::vector<float>* goodclusters_E = nullptr;
0224   std::vector<float>* goodclusters_Ecore = nullptr;
0225   std::vector<float>* goodclusters_Eta = nullptr;
0226   std::vector<float>* goodclusters_Phi = nullptr;
0227   std::vector<float>* goodclusters_pT = nullptr;
0228   std::vector<float>* goodclusters_xF = nullptr;
0229   std::vector<float>* goodclusters_Chi2 = nullptr;
0230
0231   // diphotons and cuts
0232   float min_diphotonPt = 1.0;
0233   float max_asym = 1.0;
0234   float min_deltaR = 0.0;
0235   float max_deltaR = 999.9;
0236   std::vector<float>* diphoton_E = nullptr;
0237   std::vector<float>* diphoton_M = nullptr;
0238   std::vector<float>* diphoton_Eta = nullptr;
0239   std::vector<float>* diphoton_Phi = nullptr;
0240   std::vector<float>* diphoton_pT = nullptr;
0241   std::vector<float>* diphoton_xF = nullptr;
0242   std::vector<int>* diphoton_clus1index = nullptr;
0243   std::vector<int>* diphoton_clus2index = nullptr;
0244   std::vector<float>* diphoton_deltaR = nullptr;
0245   std::vector<float>* diphoton_asym = nullptr;
0246
0247   std::pair<float, float> pi0MassRange{0.117, 0.167};
0248   std::pair<float, float> pi0BkgrLowRange{0.047, 0.097};
0249   std::pair<float, float> pi0BkgrHighRange{0.177, 0.227};
0250   std::pair<float, float> etaMassRange{0.494, 0.634};
0251   std::pair<float, float> etaBkgrLowRange{0.300, 0.400};
0252   std::pair<float, float> etaBkgrHighRange{0.700, 0.800};
0253   std::vector<Diphoton>* pi0s = nullptr;
0254   std::vector<Diphoton>* etas = nullptr;
0255   std::vector<Diphoton>* pi0Bkgr = nullptr;
0256   std::vector<Diphoton>* etaBkgr = nullptr;
0257
0258   // phi histograms for asymmetries
0259   const float PI = 3.141592;
0260   int nBins_pT = 8;
0261   float bhs_max_pT = 12.0;
0262   int nBins_xF = 8;
0263   float bhs_max_xF = 0.15;
0264   int nHistBins_phi = 16;
0265   PhiHists* pi0Hists = nullptr;
0266   PhiHists* etaHists = nullptr;
0267   PhiHists* pi0BkgrHists = nullptr;
0268   PhiHists* etaBkgrHists = nullptr;
0269   PhiHists* pi0Hists_lowEta = nullptr;
0270   PhiHists* etaHists_lowEta = nullptr;
0271   PhiHists* pi0BkgrHists_lowEta = nullptr;
0272   PhiHists* etaBkgrHists_lowEta = nullptr;
0273   PhiHists* pi0Hists_highEta = nullptr;
0274   PhiHists* etaHists_highEta = nullptr;
0275   PhiHists* pi0BkgrHists_highEta = nullptr;
0276   PhiHists* etaBkgrHists_highEta = nullptr;
0277
0278   // counters for events
0279   long int n_events_total = 0;
0280   long int mbdcoinc_withoutNandS = 0;
0281   long int n_events_mbdtrigger = 0;
0282   long int n_events_mbdtrigger_vtx1 = 0;
0283   long int n_events_mbdtrigger_vtx2 = 0;
0284   long int n_events_mbdtrigger_vtx3 = 0;
0285   long int n_events_photontrigger = 0;
0286   long int n_events_mbdvtx_first1k = 0;
0287   long int n_events_globalvtx_first1k = 0;
0288   long int first_mbdvtx = 0;
0289   long int first_globalvtx = 0;
0290   long int n_events_mbdvtx_with_mbdtrig = 0;
0291   long int n_events_mbdvtx_without_mbdtrig = 0;
0292   long int n_events_globalvtx_with_mbdtrig = 0;
0293   long int n_events_globalvtx_without_mbdtrig = 0;
0294   long int n_events_globalvtx_with_mbdvtx = 0;
0295   long int n_events_globalvtx_without_mbdvtx = 0;
0296   long int n_events_with_mbdvertex = 0;
0297   long int n_events_with_globalvertex = 0;
0298   long int n_events_with_vtx10 = 0;
0299   long int n_events_with_vtx30 = 0;
0300   long int n_events_with_vtx50 = 0;
0301   long int n_events_with_good_vertex = 0;
0302   long int n_events_positiveCaloE = 0;
0303
0304 };
0305
0306 #endif // NEUTRALMESONTSSA_H