neutralMesonTSSA/src/neutralMesonTSSA.cc

0001 //our base code
0002 #include "neutralMesonTSSA.h"
0003 #include "binnedhistset/BinnedHistSet.h"
0004
0005 //Fun4all stuff
0006 #include <fun4all/Fun4AllReturnCodes.h>
0007 #include <fun4all/Fun4AllServer.h>
0008 #include <fun4all/Fun4AllHistoManager.h>
0009 #include <phool/PHCompositeNode.h>
0010 #include <phool/getClass.h>
0011 #include <phool/phool.h>
0012 /* #include <ffaobjects/EventHeader.h> */
0013 #include <ffaobjects/RunHeader.h>
0014 #include <ffarawobjects/Gl1Packet.h>
0015
0016 //ROOT stuff
0017 #include <TH1F.h>
0018 #include <TH2F.h>
0019 #include <TFile.h>
0020 #include <TTree.h>
0021 #include <TLorentzVector.h>
0022 #include <TString.h>
0023
0024 //for emc clusters
0025 #include <calobase/RawCluster.h>
0026 #include <calobase/RawClusterContainer.h>
0027 #include <calobase/RawClusterUtility.h>
0028 #include <calobase/RawTowerGeomContainer.h>
0029 #include <calobase/RawTower.h>
0030 #include <calobase/RawTowerContainer.h>
0031 #include <calobase/TowerInfo.h>
0032 #include <calobase/TowerInfoContainer.h>
0033
0034 //spin database stuff
0035 #include <uspin/SpinDBContent.h>
0036 #include <uspin/SpinDBOutput.h>
0037
0038 // Minimum Bias
0039 /* #include <calotrigger/MinimumBiasInfo.h> */
0040
0041 //for vertex information
0042 #include <globalvertex/GlobalVertex.h>
0043 #include <globalvertex/GlobalVertexMap.h>
0044 #include <globalvertex/MbdVertex.h>
0045 #include <globalvertex/MbdVertexMap.h>
0046
0047 //truth information
0048 #include <g4main/PHG4TruthInfoContainer.h>
0049 #include <g4main/PHG4VtxPoint.h>
0050
0051 neutralMesonTSSA::neutralMesonTSSA(const std::string &name, std::string histname, std::string treename, bool isMC):
0052  SubsysReco(name),
0053  isMonteCarlo(isMC),
0054  outfilename_hists(histname),
0055  outfilename_trees(treename)
0056 {
0057   std::cout << "neutralMesonTSSA::neutralMesonTSSA(const std::string &name) Calling ctor" << std::endl;
0058 }
0059
0060 //____________________________________________________________________________..
0061 neutralMesonTSSA::~neutralMesonTSSA()
0062 {
0063   std::cout << "neutralMesonTSSA::~neutralMesonTSSA() Calling dtor" << std::endl;
0064 }
0065
0066 //____________________________________________________________________________..
0067 int neutralMesonTSSA::Init(PHCompositeNode *topNode)
0068 {
0069   std::cout << "neutralMesonTSSA::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0070
0071   outfile_hists = new TFile(outfilename_hists.c_str(), "RECREATE");
0072   outfile_hists->cd();
0073   MakeAllHists();
0074
0075   MakeVectors();
0076
0077   outfile_trees = new TFile(outfilename_trees.c_str(), "RECREATE");
0078   outfile_trees->cd();
0079
0080   tree_event = new TTree("Event", "Tree for event information");
0081   tree_event->Branch("crossingAngle", &crossingAngleIntended);
0082   tree_event->Branch("vtxz", &vtxz);
0083   tree_event->Branch("minbiastrig_live", &mbdtrigger_live);
0084   tree_event->Branch("photontrig_live", &photontrigger_live);
0085   tree_event->Branch("minbiastrig_scaled", &mbdtrigger_scaled);
0086   tree_event->Branch("photontrig_scaled", &photontrigger_scaled);
0087   tree_event->Branch("bspin", &bspin);
0088   tree_event->Branch("yspin", &yspin);
0089
0090   tree_clusters = new TTree("Clusters", "Tree for cluster information");
0091   tree_clusters->Branch("clusterE", goodclusters_E);
0092   tree_clusters->Branch("clusterEcore", goodclusters_Ecore);
0093   tree_clusters->Branch("clusterEta", goodclusters_Eta);
0094   tree_clusters->Branch("clusterPhi", goodclusters_Phi);
0095   tree_clusters->Branch("clusterpT", goodclusters_pT);
0096   tree_clusters->Branch("clusterxF", goodclusters_xF);
0097   tree_clusters->Branch("clusterChi2", goodclusters_Chi2);
0098   tree_clusters->Branch("nAllClusters", &nAllClusters);
0099
0100   tree_diphotons = new TTree("Diphotons", "Tree for diphoton information");
0101   tree_diphotons->Branch("diphotonE", diphoton_E);
0102   tree_diphotons->Branch("diphotonM", diphoton_M);
0103   tree_diphotons->Branch("diphotonEta", diphoton_Eta);
0104   tree_diphotons->Branch("diphotonPhi", diphoton_Phi);
0105   tree_diphotons->Branch("diphotonpT", diphoton_pT);
0106   tree_diphotons->Branch("diphotonxF", diphoton_xF);
0107   tree_diphotons->Branch("diphotonClusterIndex1", diphoton_clus1index);
0108   tree_diphotons->Branch("diphotonClusterIndex2", diphoton_clus2index);
0109   tree_diphotons->Branch("diphotonDeltaR", diphoton_deltaR);
0110   tree_diphotons->Branch("diphotonAsym", diphoton_asym);
0111
0112   return Fun4AllReturnCodes::EVENT_OK;
0113 }
0114
0115 //____________________________________________________________________________..
0116 int neutralMesonTSSA::InitRun(PHCompositeNode *topNode)
0117 {
0118   std::cout << "neutralMesonTSSA::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0119
0120   // Run header
0121   runHeader = findNode::getClass<RunHeader>(topNode, "RunHeader");
0122   if (!runHeader)
0123   {
0124     std::cout << PHWHERE << ":: RunHeader node missing! Skipping run XXX" << std::endl;
0125     return Fun4AllReturnCodes::ABORTRUN;
0126   }
0127
0128   GetRunNum();
0129   int bad_spin = GetSpinInfo();
0130   if (bad_spin)
0131   {
0132     return Fun4AllReturnCodes::ABORTRUN;
0133   }
0134
0135   return Fun4AllReturnCodes::EVENT_OK;
0136 }
0137
0138 //____________________________________________________________________________..
0139 int neutralMesonTSSA::process_event(PHCompositeNode *topNode)
0140 {
0141   /* std::cout << "neutralMesonTSSA::process_event(PHCompositeNode *topNode) Processing Event" << std::endl; */
0142   n_events_total++;
0143   h_nEvents->Fill(1);
0144   if (n_events_total%10000 == 0) std::cout << "Event " << n_events_total << std::endl;
0145   /* if (n_events_total < 1000) return Fun4AllReturnCodes::ABORTEVENT; */
0146   /* std::cout << "Greg info: starting process_event. n_events_total = " << n_events_total << std::endl; */
0147
0148   // First populate all the data containers
0149
0150   // GL1
0151   gl1Packet = findNode::getClass<Gl1Packet>(topNode, "GL1Packet");
0152   if (!gl1Packet)
0153   {
0154     std::cout << PHWHERE << ":: GL1Packet node missing! Skipping run " << runNum << std::endl;
0155     return Fun4AllReturnCodes::ABORTRUN;
0156   }
0157
0158   // Check for MBDNS coincidence trigger
0159   GetTrigger();
0160   if (mbdtrigger_live) h_nEvents->Fill(2);
0161   if (photontrigger_live) h_nEvents->Fill(3);
0162   if (mbdtrigger_scaled) h_nEvents->Fill(4);
0163   if (photontrigger_scaled) h_nEvents->Fill(5);
0164   /* if (mbdtrigger_live && photontrigger_live) h_nEvents->Fill(4); */
0165   if (!mbdtrigger_live) return Fun4AllReturnCodes::ABORTEVENT;
0166   /* if (!photontrigger) return Fun4AllReturnCodes::ABORTEVENT; */
0167
0168   // Information on clusters
0169   // Name of node is different in MC and RD
0170   if (isMonteCarlo) {
0171       m_clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_POS_COR_CEMC");
0172   }
0173   else {
0174       m_clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTERINFO_CEMC");
0175   }
0176   if(!m_clusterContainer)
0177   {
0178       if (isMonteCarlo) std::cout << PHWHERE << "neutralMesonTSSA::process_event - Fatal Error - CLUSTER_POS_COR_CEMC node is missing. " << std::endl;
0179       else std::cout << PHWHERE << "neutralMesonTSSA::process_event - Fatal Error - CLUSTERINFO_CEMC node is missing. " << std::endl;
0180       return Fun4AllReturnCodes::ABORTEVENT;
0181   }
0182
0183   // Raw tower information -- used for checking if total calo E > 0
0184   RawTowerContainer *towerContainer = nullptr;
0185   TowerInfoContainer *towerInfoContainer = nullptr;
0186   // again, node has different names in MC and RD
0187   if (isMonteCarlo) {
0188       towerContainer = findNode::getClass<RawTowerContainer>(topNode,"TOWER_CALIB_CEMC");
0189       if(!towerContainer) {
0190       std::cout << PHWHERE << "neutralMesonTSSA::process_event Could not find node TOWER_CALIB_CEMC"  << std::endl;
0191       return Fun4AllReturnCodes::ABORTEVENT;
0192       }
0193   }
0194   else {
0195       towerInfoContainer = findNode::getClass<TowerInfoContainer>(topNode,"TOWERINFO_CALIB_CEMC");
0196       if (!towerInfoContainer) {
0197       std::cout << PHWHERE << "neutralMesonTSSA::process_event Could not find node TOWERINFO_CALIB_CEMC"  << std::endl;
0198       return Fun4AllReturnCodes::ABORTEVENT;
0199       }
0200   }
0201
0202   // Truth information
0203   /* PHG4TruthInfoContainer *truthinfo = nullptr; */
0204   if (isMonteCarlo) {
0205       m_truthInfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0206       if(!m_truthInfo) {
0207       std::cout << PHWHERE << "neutralMesonTSSA::process_event Could not find node G4TruthInfo"  << std::endl;
0208       return Fun4AllReturnCodes::ABORTEVENT;
0209       }
0210   }
0211
0212   //Vertex information
0213   MbdVertexMap *MBDvtxContainer = findNode::getClass<MbdVertexMap>(topNode,"MbdVertexMap");
0214   if (MBDvtxContainer) {
0215       if (!MBDvtxContainer->empty()) {
0216       MbdVertex *MBDVertex= MBDvtxContainer->begin()->second;
0217       if (MBDVertex) {
0218           mbdvertex = true;
0219           n_events_with_mbdvertex++;
0220           if (mbdtrigger_live) n_events_mbdvtx_with_mbdtrig++;
0221           if (!mbdtrigger_live) n_events_mbdvtx_without_mbdtrig++;
0222           if (first_mbdvtx == 0) first_mbdvtx = n_events_total - 1;
0223           if (n_events_total < 1000) n_events_mbdvtx_first1k++;
0224       }
0225       }
0226   }
0227
0228   // Problem is MC has a PrimaryVtx but no GlobalVertex, while RD has the opposite
0229   if (isMonteCarlo) {
0230       PHG4TruthInfoContainer::VtxRange vtx_range = m_truthInfo->GetPrimaryVtxRange();
0231       PHG4TruthInfoContainer::ConstVtxIterator vtxIter = vtx_range.first;
0232       globalvertex = true;
0233       mcVtx = vtxIter->second;
0234       n_events_with_globalvertex++;
0235       vtxz = mcVtx->get_z();
0236       h_vtxz->Fill(vtxz);
0237       /* if (abs(mcVtx->get_z()) > 30.0) { */
0238       /* /1* std::cout << "Greg info: vertex z is " << mcVtx->get_z() << "\n"; *1/ */
0239       /* return Fun4AllReturnCodes::ABORTEVENT; */
0240       /* } */
0241       n_events_with_good_vertex++;
0242   }
0243   else {
0244       GlobalVertexMap *vtxContainer = findNode::getClass<GlobalVertexMap>(topNode,"GlobalVertexMap");
0245       if (!vtxContainer)
0246       {
0247       std::cout << PHWHERE << "neutralMesonTSSA::process_event - Fatal Error - GlobalVertexMap node is missing. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << std::endl;
0248       assert(vtxContainer);  // force quit
0249       return 0;
0250       }
0251       /* std::cout << "Global vertex map has size " << vtxContainer->size() << std::endl; */
0252       if (vtxContainer->empty())
0253       {
0254       // Final version:
0255       /* std::cout << PHWHERE << "neutralMesonTSSA::process_event - Fatal Error - GlobalVertexMap node is empty. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << std::endl; */
0256       return Fun4AllReturnCodes::ABORTEVENT;
0257
0258       // For testing
0259       // Use (0,0,0) as the vertex
0260       /* n_events_with_globalvertex++; */
0261       /* n_events_with_good_vertex++; */
0262       }
0263
0264       //More vertex information
0265       else {
0266       gVtx = vtxContainer->begin()->second;
0267       if (!gVtx)
0268       {
0269           /* std::cout << PHWHERE << "neutralMesonTSSA::process_event Could not find vtx from vtxContainer"  << std::endl; */
0270           return Fun4AllReturnCodes::ABORTEVENT;
0271       }
0272       globalvertex = true;
0273       n_events_with_globalvertex++;
0274       vtxz = gVtx->get_z();
0275       h_vtxz->Fill(vtxz);
0276       if (mbdtrigger_live) n_events_globalvtx_with_mbdtrig++;
0277       if (!mbdtrigger_live) n_events_globalvtx_without_mbdtrig++;
0278       if (mbdvertex) n_events_globalvtx_with_mbdvtx++;
0279       if (!mbdvertex) n_events_globalvtx_without_mbdvtx++;
0280       if (first_globalvtx == 0) first_globalvtx = n_events_total - 1;
0281       if (n_events_total < 1000) n_events_globalvtx_first1k++;
0282       // Require vertex to be within 10cm of 0
0283       if (abs(gVtx->get_z()) > 30.0) {
0284           /* std::cout << PHWHERE << ":: Vertex |z| > 30cm, skipping event" << std::endl; */
0285           /* return Fun4AllReturnCodes::ABORTEVENT; */
0286       }
0287       if (abs(gVtx->get_z()) < 10.0) {
0288           n_events_with_vtx10++;
0289       }
0290       if (abs(gVtx->get_z()) < 30.0) {
0291           n_events_with_vtx30++;
0292           n_events_with_good_vertex++;
0293       }
0294       if (abs(gVtx->get_z()) < 50.0) {
0295           n_events_with_vtx50++;
0296       }
0297       }
0298   }
0299   if (globalvertex) h_nEvents->Fill(6);
0300
0301   // Check that total calo E > 0
0302   if (isMonteCarlo) {
0303       RawTowerContainer::ConstRange tower_range = towerContainer -> getTowers();
0304       double totalCaloE = 0;
0305       for(RawTowerContainer::ConstIterator tower_iter = tower_range.first; tower_iter!= tower_range.second; tower_iter++)
0306       {
0307       double energy = tower_iter -> second -> get_energy();
0308       totalCaloE += energy;
0309       }
0310       if (totalCaloE < 0) {
0311       /* std::cout << PHWHERE << ":: Total EMCal energy < 0, skipping event" << std::endl; */
0312       return Fun4AllReturnCodes::ABORTEVENT;
0313       }
0314       n_events_positiveCaloE++;
0315   }
0316   else {
0317       unsigned int tower_range = towerInfoContainer->size();
0318       double totalCaloE = 0;
0319       for(unsigned int iter = 0; iter < tower_range; iter++)
0320       {
0321       TowerInfo* tower = towerInfoContainer->get_tower_at_channel(iter);
0322       // check if tower is good
0323       if(!tower->get_isGood()) continue;
0324       double energy = tower->get_energy();
0325       totalCaloE += energy;
0326       // eta,phi occupancy
0327       unsigned int towerkey = towerInfoContainer->encode_key(iter);
0328       int ieta = towerInfoContainer->getTowerEtaBin(towerkey);
0329       int iphi = towerInfoContainer->getTowerPhiBin(towerkey);
0330       if (energy > 0.5) h_towerEta_Phi_500MeV->Fill(ieta, iphi);
0331       if (energy > 0.8) h_towerEta_Phi_800MeV->Fill(ieta, iphi);
0332       }
0333       if (totalCaloE < 0) {
0334       /* std::cout << PHWHERE << ":: Total EMCal energy < 0, skipping event" << std::endl; */
0335       return Fun4AllReturnCodes::ABORTEVENT;
0336       }
0337       n_events_positiveCaloE++;
0338   }
0339
0340   // Event is passes all checks. Now call the functions to do the analysis
0341   /* std::cout << "Greg info: Getting bunch number" << std::endl; */
0342   GetBunchNum();
0343   /* std::cout << "Greg info: Getting spins" << std::endl; */
0344   GetSpins();
0345   /* std::cout << "Greg info: Getting good clusters. n_events_total = " << n_events_total << std::endl; */
0346   FindGoodClusters();
0347   /* std::cout << "Greg info: Getting diphotons" << std::endl; */
0348   FindDiphotons();
0349   /* std::cout << "Greg info: Filling phi hists" << std::endl; */
0350   /* FillAllPhiHists(); */
0351
0352   // Fill the trees
0353   tree_event->Fill();
0354   tree_clusters->Fill();
0355   tree_diphotons->Fill();
0356
0357   return Fun4AllReturnCodes::EVENT_OK;
0358 }
0359
0360 //____________________________________________________________________________..
0361 int neutralMesonTSSA::ResetEvent(PHCompositeNode *topNode)
0362 {
0363   /* std::cout << "neutralMesonTSSA::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl; */
0364
0365   // Clear the data containers, triggers and vertex info
0366   m_clusterContainer = nullptr;
0367   m_truthInfo = nullptr;
0368   gVtx = nullptr;
0369   mcVtx = nullptr;
0370   vtxz = -9999999.9;
0371   gl1Packet = nullptr;
0372   mbdNtrigger = false;
0373   mbdStrigger = false;
0374   mbdtrigger_live = false;
0375   photontrigger_live = false;
0376   mbdtrigger_scaled = false;
0377   photontrigger_scaled = false;
0378   mbdvertex = false;
0379   globalvertex = false;
0380   bspin = 0;
0381   yspin = 0;
0382   nAllClusters = 0;
0383   nGoodClusters = 0;
0384
0385   // Clear the vectors
0386   ClearVectors();
0387
0388   return Fun4AllReturnCodes::EVENT_OK;
0389 }
0390
0391 //____________________________________________________________________________..
0392 int neutralMesonTSSA::EndRun(const int runnumber)
0393 {
0394   std::cout << "neutralMesonTSSA::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0395
0396   // Clear the RunHeader
0397   runHeader = nullptr;
0398   return Fun4AllReturnCodes::EVENT_OK;
0399 }
0400
0401 //____________________________________________________________________________..
0402 int neutralMesonTSSA::End(PHCompositeNode *topNode)
0403 {
0404   std::cout << "neutralMesonTSSA::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0405   std::cout << "Processed " << n_events_total << " total events.\n";
0406   std::cout << "" << n_events_photontrigger << " events with photon > 3GeV trigger\n";
0407   std::cout << "" << n_events_mbdtrigger << " events with MBDN&S trigger\n";
0408   /* std::cout << "\t(" << mbdcoinc_withoutNandS << " *without* MBDN AND MBDS individually)\n"; */
0409   /* std::cout << "\t" << n_events_mbdtrigger_vtx1 << " with MBD trigger, |z_vtx| < T1\n"; */
0410   /* std::cout << "\t" << n_events_mbdtrigger_vtx2 << " with MBD trigger, |z_vtx| < T2\n"; */
0411   /* std::cout << "\t" << n_events_mbdtrigger_vtx3 << " events with MBD trigger, |z_vtx| < T3\n\n"; */
0412   std::cout << "" << n_events_with_mbdvertex << " events with MBDVertex (first event is " << first_mbdvtx << ")\n";
0413   std::cout << "\t" << n_events_mbdvtx_with_mbdtrig << " with MBDN&S trigger\n";
0414   std::cout << "\t" << n_events_mbdvtx_without_mbdtrig << " *without* MBDN&S trigger\n";
0415   std::cout << "" << n_events_with_globalvertex << " events with GlobalVertex (first event is " << first_globalvtx << ")\n";
0416   std::cout << "\t" << n_events_globalvtx_with_mbdvtx << " with MbdVertex\n";
0417   std::cout << "\t" << n_events_globalvtx_without_mbdvtx << " *without* MbdVertex\n";
0418   std::cout << "\t" << n_events_globalvtx_with_mbdtrig << " with MBDN&S trigger\n";
0419   std::cout << "\t" << n_events_globalvtx_without_mbdtrig << " *without* MBDN&S trigger\n";
0420   /* std::cout << "" << n_events_with_good_vertex << " events with |z_vtx| < 30cm.\n"; */
0421   std::cout << "" << n_events_with_vtx10 << " events with GlobalVertex |z| < 10cm.\n";
0422   std::cout << "" << n_events_with_vtx30 << " events with GlobalVertex |z| < 30cm.\n";
0423   std::cout << "" << n_events_with_vtx50 << " events with GlobalVertex |z| < 50cm.\n";
0424   std::cout << "" << n_events_positiveCaloE << " events with positive calo E.\n";
0425
0426   outfile_hists->Write();
0427   outfile_hists->Close();
0428   outfile_trees->Write();
0429   outfile_trees->Close();
0430
0431   DeleteStuff();
0432
0433   return Fun4AllReturnCodes::EVENT_OK;
0434 }
0435
0436 //____________________________________________________________________________..
0437 int neutralMesonTSSA::Reset(PHCompositeNode *topNode)
0438 {
0439  std::cout << "neutralMesonTSSA::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0440   return Fun4AllReturnCodes::EVENT_OK;
0441 }
0442
0443 //____________________________________________________________________________..
0444 void neutralMesonTSSA::Print(const std::string &what) const
0445 {
0446   std::cout << "neutralMesonTSSA::Print(const std::string &what) const Printing info for " << what << std::endl;
0447 }
0448
0449 float neutralMesonTSSA::get_min_clusterE()
0450 {
0451     return min_clusterE;
0452 }
0453
0454 void neutralMesonTSSA::set_min_clusterE(float Emin)
0455 {
0456     min_clusterE = Emin;
0457 }
0458
0459 float neutralMesonTSSA::get_max_clusterChi2()
0460 {
0461     return max_clusterChi2;
0462 }
0463
0464 void neutralMesonTSSA::set_max_clusterChi2(float Chi2max)
0465 {
0466     max_clusterChi2 = Chi2max;
0467 }
0468
0469 bool neutralMesonTSSA::InRange(float mass, std::pair<float, float> range)
0470 {
0471     bool ret = false;
0472     if ( (mass > range.first) && (mass < range.second) ) ret = true;
0473     return ret;
0474 }
0475
0476 void neutralMesonTSSA::MakePhiHists(std::string which)
0477 {
0478     PhiHists* hists = new PhiHists();
0479     std::string nameprefix, titlewhich, titlebeam;
0480     nameprefix = "h_" + which + "_";
0481
0482     if (which == "pi0") {
0483     titlewhich = "#pi^{0}";
0484     pi0Hists = hists;
0485     }
0486     else if (which == "eta") {
0487     titlewhich = "#eta";
0488     etaHists = hists;
0489     }
0490     else if (which == "pi0bkgr") {
0491     titlewhich = "#pi^{0} Background";
0492     pi0BkgrHists = hists;
0493     }
0494     else if (which == "etabkgr") {
0495     titlewhich = "#eta Background";
0496     etaBkgrHists = hists;
0497     }
0498     else if (which == "pi0_lowEta") {
0499     titlewhich = "#pi^{0}, |#eta| < 0.35";
0500     pi0Hists_lowEta = hists;
0501     }
0502     else if (which == "eta_lowEta") {
0503     titlewhich = "#eta, |#eta| < 0.35";
0504     etaHists_lowEta = hists;
0505     }
0506     else if (which == "pi0bkgr_lowEta") {
0507     titlewhich = "#pi^{0} Background, |#eta| < 0.35";
0508     pi0BkgrHists_lowEta = hists;
0509     }
0510     else if (which == "etabkgr_lowEta") {
0511     titlewhich = "#eta Background, |#eta| < 0.35";
0512     etaBkgrHists_lowEta = hists;
0513     }
0514     else if (which == "pi0_highEta") {
0515     titlewhich = "#pi^{0}, |#eta| > 0.35";
0516     pi0Hists_highEta = hists;
0517     }
0518     else if (which == "eta_highEta") {
0519     titlewhich = "#eta, |#eta| > 0.35";
0520     etaHists_highEta = hists;
0521     }
0522     else if (which == "pi0bkgr_highEta") {
0523     titlewhich = "#pi^{0} Background, |#eta| > 0.35";
0524     pi0BkgrHists_highEta = hists;
0525     }
0526     else if (which == "etabkgr_highEta") {
0527     titlewhich = "#eta Background, |#eta| > 0.35";
0528     etaBkgrHists_highEta = hists;
0529     }
0530     else {
0531     std::cout << PHWHERE << ":: Invalid arguments!" << std::endl;
0532     return;
0533     }
0534
0535     std::vector<double> pTbins = {2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 10.0, 20.0};
0536     /* std::vector<double> pTbins; */
0537     /* for (int i=0; i<nBins_pT; i++) { */
0538     /* pTbins.push_back(i*(bhs_max_pT/nBins_pT)); */
0539     /* } */
0540     /* pTbins.push_back(bhs_max_pT); */
0541     std::vector<double> xFbins = {-0.15, -0.10, -0.06, -0.03, 0.0, 0.03, 0.06, 0.10, 0.15};
0542     /* std::vector<double> xFbins; */
0543     /* for (int i=0; i<nBins_xF; i++) { */
0544     /* xFbins.push_back(2*bhs_max_xF*i/nBins_xF - bhs_max_xF); */
0545     /* } */
0546     /* xFbins.push_back(bhs_max_xF); */
0547     std::vector<double> etabins = {-2.0, -1.15, -0.35, 0.0, 0.35, 1.15, 2.0};
0548     std::vector<double> vtxzbins = {-100.0, -50.0, -30.0, 0.0, 30.0, 50.0, 100.0};
0549
0550     hists->phi_pT = new BinnedHistSet(Form("%sphi_pT", nameprefix.c_str()), Form("%s #phi Distribution;#phi (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "p_{T} (GeV)", pTbins);
0551     hists->phi_pT_blue_up = new BinnedHistSet(Form("%sphi_pT_blue_up", nameprefix.c_str()), Form("%s Blue Beam Spin-Up #phi^{B} Distribution;#phi^{B} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "p_{T} (GeV)", pTbins);
0552     hists->phi_pT_blue_down = new BinnedHistSet(Form("%sphi_pT_blue_down", nameprefix.c_str()), Form("%s Blue Beam Spin-Down #phi^{B} Distribution;#phi^{B} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "p_{T} (GeV)", pTbins);
0553     hists->phi_pT_yellow_up = new BinnedHistSet(Form("%sphi_pT_yellow_up", nameprefix.c_str()), Form("%s Yellow Beam Spin-Up #phi^{Y} Distribution;#phi^{Y} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "p_{T} (GeV)", pTbins);
0554     hists->phi_pT_yellow_down = new BinnedHistSet(Form("%sphi_pT_yellow_down", nameprefix.c_str()), Form("%s Yellow Beam Spin-Down #phi^{Y} Distribution;#phi^{Y} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "p_{T} (GeV)", pTbins);
0555
0556     hists->phi_xF = new BinnedHistSet(Form("%sphi_xF", nameprefix.c_str()), Form("%s #phi Distribution;#phi (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "x_{F}", xFbins);
0557     hists->phi_xF_blue_up = new BinnedHistSet(Form("%sphi_xF_blue_up", nameprefix.c_str()), Form("%s Blue Beam Spin-Up #phi^{B} Distribution;#phi^{B} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "x_{F}", xFbins);
0558     hists->phi_xF_blue_down = new BinnedHistSet(Form("%sphi_xF_blue_down", nameprefix.c_str()), Form("%s Blue Beam Spin-Down #phi^{B} Distribution;#phi^{B} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "x_{F}", xFbins);
0559     hists->phi_xF_yellow_up = new BinnedHistSet(Form("%sphi_xF_yellow_up", nameprefix.c_str()), Form("%s Yellow Beam Spin-Up #phi^{Y} Distribution;#phi^{Y} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "x_{F}", xFbins);
0560     hists->phi_xF_yellow_down = new BinnedHistSet(Form("%sphi_xF_yellow_down", nameprefix.c_str()), Form("%s Yellow Beam Spin-Down #phi^{Y} Distribution;#phi^{Y} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "x_{F}", xFbins);
0561
0562     hists->phi_eta = new BinnedHistSet(Form("%sphi_eta", nameprefix.c_str()), Form("%s #phi Distribution;#phi (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "#eta", etabins);
0563     hists->phi_eta_blue_up = new BinnedHistSet(Form("%sphi_eta_blue_up", nameprefix.c_str()), Form("%s Blue Beam Spin-Up #phi^{B} Distribution;#phi^{B} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "#eta", etabins);
0564     hists->phi_eta_blue_down = new BinnedHistSet(Form("%sphi_eta_blue_down", nameprefix.c_str()), Form("%s Blue Beam Spin-Down #phi^{B} Distribution;#phi^{B} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "#eta", etabins);
0565     hists->phi_eta_yellow_up = new BinnedHistSet(Form("%sphi_eta_yellow_up", nameprefix.c_str()), Form("%s Yellow Beam Spin-Up #phi^{Y} Distribution;#phi^{Y} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "#eta", etabins);
0566     hists->phi_eta_yellow_down = new BinnedHistSet(Form("%sphi_eta_yellow_down", nameprefix.c_str()), Form("%s Yellow Beam Spin-Down #phi^{Y} Distribution;#phi^{Y} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "#eta", etabins);
0567
0568     hists->phi_vtxz = new BinnedHistSet(Form("%sphi_vtxz", nameprefix.c_str()), Form("%s #phi Distribution;#phi (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "z_{vtx} (cm)", vtxzbins);
0569     hists->phi_vtxz_blue_up = new BinnedHistSet(Form("%sphi_vtxz_blue_up", nameprefix.c_str()), Form("%s Blue Beam Spin-Up #phi^{B} Distribution;#phi^{B} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "z_{vtx} (cm)", vtxzbins);
0570     hists->phi_vtxz_blue_down = new BinnedHistSet(Form("%sphi_vtxz_blue_down", nameprefix.c_str()), Form("%s Blue Beam Spin-Down #phi^{B} Distribution;#phi^{B} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "z_{vtx} (cm)", vtxzbins);
0571     hists->phi_vtxz_yellow_up = new BinnedHistSet(Form("%sphi_vtxz_yellow_up", nameprefix.c_str()), Form("%s Yellow Beam Spin-Up #phi^{Y} Distribution;#phi^{Y} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "z_{vtx} (cm)", vtxzbins);
0572     hists->phi_vtxz_yellow_down = new BinnedHistSet(Form("%sphi_vtxz_yellow_down", nameprefix.c_str()), Form("%s Yellow Beam Spin-Down #phi^{Y} Distribution;#phi^{Y} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "z_{vtx} (cm)", vtxzbins);
0573
0574     hists->phi_pT->MakeHists();
0575     hists->phi_pT_blue_up->MakeHists();
0576     hists->phi_pT_blue_down->MakeHists();
0577     hists->phi_pT_yellow_up->MakeHists();
0578     hists->phi_pT_yellow_down->MakeHists();
0579
0580     hists->phi_xF->MakeHists();
0581     hists->phi_xF_blue_up->MakeHists();
0582     hists->phi_xF_blue_down->MakeHists();
0583     hists->phi_xF_yellow_up->MakeHists();
0584     hists->phi_xF_yellow_down->MakeHists();
0585
0586     hists->phi_eta->MakeHists();
0587     hists->phi_eta_blue_up->MakeHists();
0588     hists->phi_eta_blue_down->MakeHists();
0589     hists->phi_eta_yellow_up->MakeHists();
0590     hists->phi_eta_yellow_down->MakeHists();
0591
0592     hists->phi_vtxz->MakeHists();
0593     hists->phi_vtxz_blue_up->MakeHists();
0594     hists->phi_vtxz_blue_down->MakeHists();
0595     hists->phi_vtxz_yellow_up->MakeHists();
0596     hists->phi_vtxz_yellow_down->MakeHists();
0597 }
0598
0599 void neutralMesonTSSA::MakeAllHists()
0600 {
0601     // trigger & vertex
0602     // h_nEvents: 1 = all events; 2 = minbias trig; 3 = photon trig; 4 = both trig; 5 = GlobalVertex
0603     h_nEvents = new TH1F("h_nEvents", "Number of Events Processed", 6, 0.5, 6.5);
0604     h_nEvents->GetXaxis()->SetBinLabel(1, "All Events");
0605     h_nEvents->GetXaxis()->SetBinLabel(2, "Minbias Trig Live");
0606     h_nEvents->GetXaxis()->SetBinLabel(3, "Photon Trig Live");
0607     h_nEvents->GetXaxis()->SetBinLabel(4, "Minbias Trig Scaled");
0608     h_nEvents->GetXaxis()->SetBinLabel(5, "Photon Trig Scaled");
0609     h_nEvents->GetXaxis()->SetBinLabel(6, "GlobalVertex");
0610
0611     int ntowers_eta = 96;
0612     int ntowers_phi = 256;
0613     h_towerEta_Phi_500MeV = new TH2F("h_towerEta_Phi_500MeV", "Towers with E > 0.5GeV;i_{#eta};i_{#phi}", ntowers_eta, 0, ntowers_eta, ntowers_phi, 0, ntowers_phi);
0614     h_towerEta_Phi_800MeV = new TH2F("h_towerEta_Phi_800MeV", "Towers with E > 0.8GeV;i_{#eta};i_{#phi}", ntowers_eta, 0, ntowers_eta, ntowers_phi, 0, ntowers_phi);
0615
0616     int nbins_vtxz = 100;
0617     double vtxz_upper = 200.0;
0618     h_vtxz = new TH1F("h_vtxz", "Vertex z Distribution;z_{vtx} (cm);Counts", nbins_vtxz, -vtxz_upper, vtxz_upper);
0619     h_crossingAngle = new TH1F("h_crossingAngle", "Crossing Angle;Crossing Angle (mrad);# Runs", 4, 0.0, 4.0);
0620     h_crossingAngle->GetXaxis()->SetBinLabel(1, "Default -999");
0621     h_crossingAngle->GetXaxis()->SetBinLabel(2, "-1.5");
0622     h_crossingAngle->GetXaxis()->SetBinLabel(3, "0");
0623     h_crossingAngle->GetXaxis()->SetBinLabel(4, "+1.5");
0624
0625     // all clusters
0626     int nbins_etaphi = 200;
0627     double eta_upper = 2.00;
0628     double phi_upper = PI;
0629     int nbins_pT = 100;
0630     int nbins_xF = 100;
0631     double pT_upper = 20.0;
0632     double xF_upper = 0.20;
0633     h_nAllClusters = new TH1F("h_nAllClusters", "Total Number of Clusters per Event;# Clusters;Counts", 500, 800, 2800);
0634     h_nGoodClusters = new TH1F("h_nGoodClusters", "Number of \"Good\" Clusters per Event;# Clusters;Counts", 10, -0.5, 9.5);
0635     h_clusterE = new TH1F("h_clusterE", "Cluster Energy Distribution;Cluster E (GeV);Counts", nbins_pT, 0.0, pT_upper);
0636     h_clusterEta = new TH1F("h_clusterEta", "Cluster #eta Distribution;Cluster #eta;Counts", nbins_etaphi, -eta_upper, eta_upper);
0637     h_clusterVtxz_Eta = new TH2F("h_clusterVtxz_Eta", "Cluster #eta v. Vertex z;Vertex z (cm);Cluster #eta", nbins_etaphi, -vtxz_upper, vtxz_upper, nbins_etaphi, -eta_upper, eta_upper);
0638     h_clusterPhi = new TH1F("h_clusterPhi", "Cluster #phi Distribution;Cluster #phi;Counts", nbins_etaphi, -phi_upper, phi_upper);
0639     h_clusterEta_Phi = new TH2F("h_clusterEta_Phi", "Cluster Position;Cluster #eta;Cluster #phi (rad)", nbins_etaphi, -eta_upper, eta_upper, nbins_etaphi, -phi_upper, phi_upper);
0640     h_clusterpT = new TH1F("h_clusterpT", "Cluster Transverse Momentum Distribution;Cluster p_{T} (GeV);Counts", nbins_pT, 0.0, pT_upper);
0641     h_clusterxF = new TH1F("h_clusterxF", "Cluster x_{F} Distribution;Cluster x_{F};Counts", nbins_xF, -1.0*xF_upper, xF_upper);
0642     h_clusterpT_xF = new TH2F("h_clusterpT_xF", "Cluster x_{F} vs p_{T};Cluster p_{T} (GeV);Cluster x_{F}", nbins_pT, 0.0, pT_upper, nbins_xF, -1.0*xF_upper, xF_upper);
0643     h_clusterChi2 = new TH1F("h_clusterChi2", "Cluster #chi^{2} Distribution;Cluster #chi^{2};Counts", 200, 0.0, 50.0);
0644     h_clusterChi2zoomed = new TH1F("h_clusterChi2zoomed", "Cluster #chi^{2} Distribution;Cluster #chi^{2};Counts", 200, 0.0, 5.0);
0645     h_mesonClusterChi2 = new TH1F("h_mesonClusterChi2", "Cluster #chi^{2} for #pi^{0} and #eta Clusters;Cluster #chi^{2};Counts", 200, 0.0, 5.0);
0646
0647     // good clusters
0648     h_goodClusterE = new TH1F("h_goodClusterE", "Energy Distribution of \"Good\" Clusters;Cluster E (GeV);Counts", nbins_pT, 0.0, pT_upper);
0649     h_goodClusterEta = new TH1F("h_goodClusterEta", "#eta Distribution of \"Good\" Clusters;Cluster #eta;Counts", nbins_etaphi, -eta_upper, eta_upper);
0650     h_goodClusterVtxz_Eta = new TH2F("h_goodClusterVtxz_Eta", "\"Good\" Cluster #eta v. Vertex z;Vertex z (cm);Cluster #eta", nbins_etaphi, -vtxz_upper, vtxz_upper, nbins_etaphi, -eta_upper, eta_upper);
0651     h_goodClusterPhi = new TH1F("h_goodClusterPhi", "#phi Distribution of \"Good\" Clusters;Cluster #phi;Counts", nbins_etaphi, -phi_upper, phi_upper);
0652     h_goodClusterEta_Phi = new TH2F("h_goodClusterEta_Phi", "Angular Position of \"Good\" Clusters;Cluster #eta;Cluster #phi (rad)", nbins_etaphi, -eta_upper, eta_upper, nbins_etaphi, -phi_upper, phi_upper);
0653     h_goodClusterpT = new TH1F("h_goodClusterpT", "Transverse Momentum Distribution of \"Good\" Clusters;Cluster p_{T} (GeV);Counts", nbins_pT, 0.0, pT_upper);
0654     h_goodClusterxF = new TH1F("h_goodClusterxF", "x_{F} Distribution of \"Good\" Clusters;Cluster x_{F};Counts", nbins_xF, -1.0*xF_upper, xF_upper);
0655     h_goodClusterpT_xF = new TH2F("h_goodClusterpT_xF", "\"Good\" Cluster x_{F} vs p_{T};Cluster p_{T} (GeV);Cluster x_{F}", nbins_pT, 0.0, pT_upper, nbins_xF, -1.0*xF_upper, xF_upper);
0656
0657     // diphotons
0658     h_nDiphotons = new TH1F("h_nDiphotons", "Number of Diphotons per Event;# #gamma #gamma pairs;Counts", 100, -0.5, 99.5);
0659     h_nRecoPi0s = new TH1F("h_nRecoPi0s", "Number of Reconstructed #pi^{0}s per Event;# #pi^{0}s;Counts", 10, -0.5, 9.5);
0660     h_nRecoEtas = new TH1F("h_nRecoEtas", "Number of Reconstructed #eta_{}s per Event;# #eta_{}s;Counts", 10, -0.5, 9.5);
0661     int nbins_mass = 100;
0662     h_diphotonMass = new TH1F("h_diphotonMass", "Diphoton Mass Distribution;Diphoton Mass (GeV);Counts", nbins_mass, 0.0, 1.0);
0663     h_diphotonE = new TH1F("h_diphotonE", "Diphoton Energy Distribution;Diphoton E (GeV);Counts", nbins_pT, 0.0, pT_upper);
0664     h_diphotonEta = new TH1F("h_diphotonEta", "Diphoton #eta Distribution;Diphoton #eta;Counts", nbins_etaphi, -eta_upper, eta_upper);
0665     h_diphotonVtxz_Eta = new TH2F("h_diphotonVtxz_Eta", "Diphoton #eta v. Vertex z;Vertex z (cm);Diphoton #eta", nbins_etaphi, -vtxz_upper, vtxz_upper, nbins_etaphi, -eta_upper, eta_upper);
0666     h_diphotonPhi = new TH1F("h_diphotonPhi", "Diphoton #phi Distribution;Diphoton #phi;Counts", nbins_etaphi, -phi_upper, phi_upper);
0667     h_diphotonEta_Phi = new TH2F("h_diphotonEta_Phi", "Diphoton Position;Diphoton #eta;Diphoton #phi (rad)", nbins_etaphi, -eta_upper, eta_upper, nbins_etaphi, -phi_upper, phi_upper);
0668     h_diphotonpT = new TH1F("h_diphotonpT", "Diphoton Transverse Momentum Distribution;Diphoton p_{T} (GeV);Counts", nbins_pT, 0.0, pT_upper);
0669     h_diphotonxF = new TH1F("h_diphotonxF", "Diphoton x_{F} Distribution;Diphoton x_{F};Counts", nbins_xF, -1.0*xF_upper, xF_upper);
0670     h_diphotonpT_xF = new TH2F("h_diphotonpT_xF", "Diphoton x_{F} vs p_{T};Diphoton p_{T} (GeV);Diphoton x_{F}", nbins_pT, 0.0, pT_upper, nbins_xF, -1.0*xF_upper, xF_upper);
0671     h_diphotonAsym = new TH1F("h_diphotonAsym", "Diphoton Pair Asymmetry Distribution;#alpha;Counts", nbins_xF, 0.0, 1.0);
0672     h_diphotonDeltaR = new TH1F("h_diphotonDeltaR", "Diphoton #Delta R Distribution;Diphoton #Delta R;Counts", nbins_xF, 0.0, PI);
0673     h_diphotonAsym_DeltaR = new TH2F("h_diphotonAsym_DeltaR", "Diphoton #Delta R vs Pair Asymmetry;#alpha;#Delta R", nbins_xF, 0.0, 1.0, nbins_xF, 0.0, PI);
0674
0675     std::vector<double> pTbins;
0676     std::vector<double> xFbins;
0677     int nbins_bhs = 20;
0678     for (int i=0; i<nbins_bhs; i++) {
0679     pTbins.push_back(i*(pT_upper/nbins_bhs));
0680     xFbins.push_back(2*xF_upper*i/nbins_bhs - xF_upper);
0681     }
0682     pTbins.push_back(pT_upper);
0683     xFbins.push_back(xF_upper);
0684
0685     bhs_diphotonMass_pT = new BinnedHistSet("h_diphotonMass_pT", "Diphoton Mass;Mass (GeV);Counts", nbins_mass, 0.0, 1.0, "p_{T} (GeV)", pTbins);
0686     bhs_diphotonMass_xF = new BinnedHistSet("h_diphotonMass_xF", "Diphoton Mass;Mass (GeV);Counts", nbins_mass, 0.0, 1.0, "x_{F}", xFbins);
0687     bhs_diphotonMass_pT->MakeHists();
0688     bhs_diphotonMass_xF->MakeHists();
0689
0690     /*
0691     MakePhiHists("pi0");
0692     MakePhiHists("eta");
0693     MakePhiHists("pi0bkgr");
0694     MakePhiHists("etabkgr");
0695     MakePhiHists("pi0_lowEta");
0696     MakePhiHists("eta_lowEta");
0697     MakePhiHists("pi0bkgr_lowEta");
0698     MakePhiHists("etabkgr_lowEta");
0699     MakePhiHists("pi0_highEta");
0700     MakePhiHists("eta_highEta");
0701     MakePhiHists("pi0bkgr_highEta");
0702     MakePhiHists("etabkgr_highEta");
0703     */
0704 }
0705
0706 void neutralMesonTSSA::MakeVectors()
0707 {
0708     goodclusters_E = new std::vector<float>;
0709     goodclusters_Ecore = new std::vector<float>;
0710     goodclusters_Eta = new std::vector<float>;
0711     goodclusters_Phi = new std::vector<float>;
0712     goodclusters_pT = new std::vector<float>;
0713     goodclusters_xF = new std::vector<float>;
0714     goodclusters_Chi2 = new std::vector<float>;
0715
0716     diphoton_E = new std::vector<float>;
0717     diphoton_M = new std::vector<float>;
0718     diphoton_Eta = new std::vector<float>;
0719     diphoton_Phi = new std::vector<float>;
0720     diphoton_pT = new std::vector<float>;
0721     diphoton_xF = new std::vector<float>;
0722     diphoton_clus1index = new std::vector<int>;
0723     diphoton_clus2index = new std::vector<int>;
0724     diphoton_deltaR = new std::vector<float>;
0725     diphoton_asym = new std::vector<float>;
0726
0727     /* pi0s = new std::vector<Diphoton>; */
0728     /* etas = new std::vector<Diphoton>; */
0729     /* pi0Bkgr = new std::vector<Diphoton>; */
0730     /* etaBkgr = new std::vector<Diphoton>; */
0731 }
0732
0733 void neutralMesonTSSA::GetRunNum()
0734 {
0735     if (!runHeader)
0736     {
0737     std::cout << PHWHERE << ":: Missing RunHeader! Exiting!" << std::endl;
0738     exit(1);
0739     }
0740     runNum = runHeader->get_RunNumber();
0741     /* std::cout << "Run number is " << runNum << std::endl; */
0742 }
0743
0744 int neutralMesonTSSA::GetSpinInfo()
0745 {
0746     int spinDB_status = 0;
0747     // Get the spin patterns from the spin DB
0748     //  0xffff is the qa_level from XingShiftCal //////
0749     unsigned int qa_level = 0xffff;
0750     SpinDBContent spin_cont;
0751     SpinDBOutput spin_out("phnxrc");
0752
0753     spin_out.StoreDBContent(runNum,runNum,qa_level);
0754     spin_out.GetDBContentStore(spin_cont,runNum);
0755
0756     // Get crossing shift
0757     crossingShift = spin_cont.GetCrossingShift();
0758     if (crossingShift == -999)
0759     {
0760     std::cout << "Warning: found crossing shift = -999 from Spin DB." << std::endl;
0761     crossingShift = 0;
0762     }
0763     std::cout << "Crossing shift: " << crossingShift << std::endl;
0764
0765     // Get spin patterns
0766     std::cout << "Blue spin pattern: [";
0767     for (int i = 0; i < NBUNCHES; i++)
0768     {
0769     spinPatternBlue[i] = spin_cont.GetSpinPatternBlue(i);
0770     std::cout << spinPatternBlue[i];
0771     if (i < 119)std::cout << ", ";
0772     }
0773     std::cout << "]" << std::endl;
0774
0775     std::cout << "Yellow spin pattern: [";
0776     for (int i = 0; i < NBUNCHES; i++)
0777     {
0778     spinPatternYellow[i] = spin_cont.GetSpinPatternYellow(i);
0779     std::cout << spinPatternYellow[i];
0780     if (i < 119)std::cout << ", ";
0781     }
0782     std::cout << "]" << std::endl;
0783
0784     // Get beam polarizations
0785     float bpolerr, ypolerr;
0786     spin_cont.GetPolarizationBlue(0, polBlue, bpolerr);
0787     spin_cont.GetPolarizationYellow(0, polYellow, ypolerr);
0788
0789     // Get GL1P scalers
0790     std::cout << "MBDNS GL1p scalers: [";
0791     for (int i = 0; i < NBUNCHES; i++)
0792     {
0793         gl1pScalers[i] = spin_cont.GetScalerMbdNoCut(i);
0794     std::cout << gl1pScalers[i];
0795     if (i < 119) std::cout << ", ";
0796     }
0797     std::cout << "]" << std::endl;
0798
0799     // Calculate luminosities
0800     for (int i = 0; i < NBUNCHES; i++)
0801     {
0802         int bsp = spinPatternBlue[i];
0803         int ysp = spinPatternYellow[i];
0804     if (bsp == 1) lumiUpBlue += gl1pScalers[i];
0805     if (bsp == -1) lumiDownBlue += gl1pScalers[i];
0806     if (ysp == 1) lumiUpYellow += gl1pScalers[i];
0807     if (ysp == -1) lumiDownYellow += gl1pScalers[i];
0808     }
0809
0810     // Print run number, relative luminosity & polarization, total luminosity
0811     std::cout << "Luminosity info: run#,bPol,bLumiUp,bLumiDown,yPol,yLumiUp,yLumiDown" << std::endl;
0812     std::cout << Form("%d,%f,%f,%f,%f,%f,%f", runNum, polBlue/100.0, lumiUpBlue, lumiDownBlue, polYellow/100.0, lumiUpYellow, lumiDownYellow) << std::endl;
0813
0814     // Get crossing angle
0815     crossingAngle = spin_cont.GetCrossAngle();
0816     // set "intended" angle to 0 or +-1.5 depending on which is closest
0817     if (crossingAngle < -0.75) {
0818     crossingAngleIntended = -1.5;
0819     h_crossingAngle->Fill(1);
0820     }
0821     else if (crossingAngle > 0.75) {
0822     crossingAngleIntended = 1.5;
0823     h_crossingAngle->Fill(3);
0824     }
0825     else {
0826     crossingAngleIntended = 0.0;
0827     h_crossingAngle->Fill(2);
0828     }
0829     if (crossingAngle == -999) {
0830     std::cout << "Warning: crossing angle set to -999. Saving 0 instead" << std::endl;
0831     crossingAngleIntended = -999;
0832     h_crossingAngle->Fill(0);
0833     }
0834
0835     // Check if spin info is valid
0836     if (spinPatternYellow[0] == -999) spinDB_status = 1;
0837     if (lumiUpBlue == 0) spinDB_status = 2;
0838     int badRunFlag = spin_cont.GetBadRunFlag();
0839     /* std::cout << "badRunFlag = " << badRunFlag << std::endl; */
0840     if (badRunFlag) spinDB_status = 3;
0841
0842     if (spinDB_status)
0843     {
0844     std::cout << PHWHERE << ":: Run number " << runNum << " has bad spin DB info! Skipping this run!" << std::endl;
0845     if (spinDB_status == 1) std::cout << "Spin pattern not stored" << std::endl;
0846     if (spinDB_status == 2) std::cout << "GL1P scalers empty" << std::endl;
0847     if (spinDB_status == 3) std::cout << "BadRunFlag set" << std::endl;
0848     /* std::cout << PHWHERE << ":: Run number " << runNum << " not found in spin DB! Exiting!" << std::endl; */
0849     /* exit(1); */
0850     }
0851     return spinDB_status;
0852 }
0853
0854 void neutralMesonTSSA::GetBunchNum()
0855 {
0856   if (gl1Packet)
0857   {
0858     bunchNum = gl1Packet->getBunchNumber();
0859     sphenixBunch = (bunchNum + crossingShift)%NBUNCHES;
0860   }
0861   else
0862   {
0863     std::cout << "GL1 missing!" << std::endl;
0864     exit(1);
0865   }
0866 }
0867
0868 void neutralMesonTSSA::GetSpins()
0869 {
0870     bspin = spinPatternBlue[sphenixBunch];
0871     yspin = spinPatternYellow[sphenixBunch];
0872 }
0873
0874 /* void neutralMesonTSSA::CountLumi() */
0875 /* { */
0876 /*     float bLumi = 0; // Replace with GL1 Scaler value */
0877 /*     if (bspin == 1) lumiUpBlue += bLumi; */
0878 /*     if (bspin == -1) lumiDownBlue += bLumi; */
0879 /*     float yLumi = 0; // Replace with GL1 Scaler value */
0880 /*     if (yspin == 1) lumiUpYellow += yLumi; */
0881 /*     if (yspin == -1) lumiDownYellow += yLumi; */
0882 /* } */
0883
0884 void neutralMesonTSSA::FindGoodClusters()
0885 {
0886     RawClusterContainer::ConstRange clusterRange = m_clusterContainer -> getClusters();
0887     RawClusterContainer::ConstIterator clusterIter;
0888     /* std::cout << "\n\nBeginning cluster loop\n"; */
0889     /* std::cout << "cluster container size is " << m_clusterContainer->size() << std::endl; */
0890     /* double vtxz = -999999.9; */
0891     /* if (gVtx) vtxz = gVtx->get_z(); */
0892
0893     for (clusterIter = clusterRange.first; clusterIter != clusterRange.second; clusterIter++)
0894     {
0895     /* std::cout << "clusterIter->first = " << clusterIter->first << std::endl; */
0896     RawCluster *recoCluster = clusterIter -> second;
0897     /* std::cout << "recoCluster = " << recoCluster << std::endl; */
0898     CLHEP::Hep3Vector vertex;
0899     if (isMonteCarlo) {
0900         /* std::cout << "Setting MC vertex, mcVtx = " << mcVtx << std::endl; */
0901         vertex = CLHEP::Hep3Vector(mcVtx->get_x(), mcVtx->get_y(), mcVtx->get_z());
0902     }
0903     else {
0904         /* std::cout << "Setting RD vertex, gVtx = " << gVtx << ", identify():" << std::endl; */
0905         if (gVtx) {
0906         /* int isValid = gVtx->isValid(); */
0907         /* std::cout << "Vertex isValid = " << isValid << std::endl; */
0908         /* float vtx_z = gVtx->get_z(); */
0909         /* std::cout << "Vertex z = " << vtx_z << std::endl; */
0910         /* gVtx->identify(); */
0911             vertex = CLHEP::Hep3Vector(gVtx->get_x(), gVtx->get_y(), gVtx->get_z());
0912         }
0913         else {
0914         /* vertex = CLHEP::Hep3Vector(0,0,0); */
0915         std::cout << "Warning: gVtx is empty!" << std::endl;
0916         continue;
0917         }
0918         /* std::cout << "Done setting RD vertex, gVtx = " << gVtx << std::endl; */
0919     }
0920     /* std::cout << "Got vertex" << std::endl; */
0921     CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
0922     /* std::cout << "E_vec_cluster" << std::endl; */
0923     float clusE = recoCluster->get_energy();
0924     float clusEcore = recoCluster->get_ecore();
0925     float clus_eta = E_vec_cluster.pseudoRapidity();
0926     float clus_phi = E_vec_cluster.phi();
0927     float clus_chi2 = recoCluster->get_chi2();
0928     float clus_pT = clusEcore/TMath::CosH(clus_eta);
0929     float clus_xF = 2.0*(clusEcore*TMath::TanH(clus_eta))/200.0;
0930
0931     nAllClusters++;
0932     h_clusterE->Fill(clusE);
0933     h_clusterEta->Fill(clus_eta);
0934     h_clusterVtxz_Eta->Fill(vtxz, clus_eta);
0935     h_clusterPhi->Fill(clus_phi);
0936     h_clusterEta_Phi->Fill(clus_eta, clus_phi);
0937     h_clusterpT->Fill(clus_pT);
0938     h_clusterxF->Fill(clus_xF);
0939     h_clusterpT_xF->Fill(clus_pT, clus_xF);
0940     h_clusterChi2->Fill(clus_chi2);
0941     h_clusterChi2zoomed->Fill(clus_chi2);
0942
0943     /* std::cout << "Applying cuts" << std::endl; */
0944     if (clusEcore < min_clusterE) continue;
0945     if (clusEcore > max_clusterE) continue;
0946     if (clus_chi2 > max_clusterChi2) continue;
0947
0948     nGoodClusters++;
0949     h_goodClusterE->Fill(clusE);
0950     h_goodClusterEta->Fill(clus_eta);
0951     h_goodClusterVtxz_Eta->Fill(vtxz, clus_eta);
0952     h_goodClusterPhi->Fill(clus_phi);
0953     h_goodClusterEta_Phi->Fill(clus_eta, clus_phi);
0954     h_goodClusterpT->Fill(clus_pT);
0955     h_goodClusterxF->Fill(clus_xF);
0956     h_goodClusterpT_xF->Fill(clus_pT, clus_xF);
0957     /* std::cout << "Populating goodcluster vectors" << std::endl; */
0958     goodclusters_E->push_back(clusE);
0959     goodclusters_Ecore->push_back(clusEcore);
0960     goodclusters_Eta->push_back(clus_eta);
0961     goodclusters_Phi->push_back(clus_phi);
0962     goodclusters_pT->push_back(clus_pT);
0963     goodclusters_xF->push_back(clus_xF);
0964     goodclusters_Chi2->push_back(clus_chi2);
0965     /* std::cout << "Done" << std::endl; */
0966     }
0967     h_nAllClusters->Fill(nAllClusters);
0968     h_nGoodClusters->Fill(nGoodClusters);
0969     return;
0970 }
0971
0972 void neutralMesonTSSA::FindDiphotons()
0973 {
0974     int nDiphotons = 0;
0975     int nRecoPi0s = 0;
0976     int nRecoEtas = 0;
0977
0978     /* double vtxz = -999999.9; */
0979     /* if (gVtx) vtxz = gVtx->get_z(); */
0980
0981     if (nGoodClusters < 2) return;
0982
0983     for (int i=0; i<(nGoodClusters-1); i++) {
0984     float E1 = goodclusters_Ecore->at(i);
0985     for (int j=(i+1); j<nGoodClusters; j++) {
0986         float E2 = goodclusters_Ecore->at(j);
0987         nDiphotons++;
0988         int lead_idx = i;
0989         int sub_idx = j;
0990         if (E1 < E2) {
0991         lead_idx = j;
0992         sub_idx = i;
0993         }
0994
0995         double leadE = goodclusters_Ecore->at(lead_idx);
0996         double leadEta = goodclusters_Eta->at(lead_idx);
0997         double leadPhi = goodclusters_Phi->at(lead_idx);
0998         double leadChi2 = goodclusters_Chi2->at(lead_idx);
0999         double subE = goodclusters_Ecore->at(sub_idx);
1000         double subEta = goodclusters_Eta->at(sub_idx);
1001         double subPhi = goodclusters_Phi->at(sub_idx);
1002         double subChi2 = goodclusters_Chi2->at(sub_idx);
1003         double asymmetry = (leadE - subE)/(leadE + subE);
1004         if (asymmetry > max_asym) continue;
1005
1006         TLorentzVector lead, sub;
1007         float leadpT = leadE/TMath::CosH(leadEta);
1008         float subpT = subE/TMath::CosH(subEta);
1009         lead.SetPtEtaPhiE(leadpT, leadEta, leadPhi, leadE);
1010         sub.SetPtEtaPhiE(subpT, subEta, subPhi, subE);
1011         TLorentzVector diphoton = lead + sub;
1012         if (diphoton.Perp() < min_diphotonPt) continue;
1013         double xF = 2.0*diphoton.Pz()/200.0;
1014         double dR = lead.DeltaR(sub);
1015         if (dR < min_deltaR) continue;
1016         if (dR > max_deltaR) continue;
1017
1018         h_diphotonE->Fill(diphoton.E());
1019         h_diphotonMass->Fill(diphoton.M());
1020         h_diphotonEta->Fill(diphoton.Eta());
1021         h_diphotonVtxz_Eta->Fill(vtxz, diphoton.Eta());
1022         h_diphotonPhi->Fill(diphoton.Phi());
1023         h_diphotonEta_Phi->Fill(diphoton.Eta(), diphoton.Phi());
1024         h_diphotonpT->Fill(diphoton.Perp());
1025         h_diphotonxF->Fill(xF);
1026         h_diphotonpT_xF->Fill(diphoton.Perp(), xF);
1027         h_diphotonAsym->Fill(asymmetry);
1028         h_diphotonDeltaR->Fill(dR);
1029         h_diphotonAsym_DeltaR->Fill(asymmetry, dR);
1030         bhs_diphotonMass_pT->FillHists(diphoton.Perp(), diphoton.M());
1031         bhs_diphotonMass_xF->FillHists(xF, diphoton.M());
1032
1033         diphoton_E->push_back(diphoton.E());
1034         diphoton_M->push_back(diphoton.M());
1035         diphoton_Eta->push_back(diphoton.Eta());
1036         diphoton_Phi->push_back(diphoton.Phi());
1037         diphoton_pT->push_back(diphoton.Perp());
1038         diphoton_xF->push_back(xF);
1039         diphoton_clus1index->push_back(lead_idx);
1040         diphoton_clus2index->push_back(sub_idx);
1041         diphoton_deltaR->push_back(dR);
1042         diphoton_asym->push_back(asymmetry);
1043
1044         /* Diphoton d; */
1045         /* d.mass = diphoton.M(); */
1046         /* d.E = diphoton.E(); */
1047         /* d.eta = diphoton.Eta(); */
1048         /* d.phi = diphoton.Phi(); */
1049         /* d.pT = diphoton.Perp(); */
1050         /* d.xF = xF; */
1051         /* d.vtxz = vtxz; */
1052
1053         /* std::vector<Diphoton>* vec = nullptr; */
1054         if (InRange(diphoton.M(), pi0MassRange)) {
1055         nRecoPi0s++;
1056         /* vec = pi0s; */
1057         h_mesonClusterChi2->Fill(leadChi2);
1058         h_mesonClusterChi2->Fill(subChi2);
1059         }
1060         if (InRange(diphoton.M(), etaMassRange)) {
1061         nRecoEtas++;
1062         /* vec = etas; */
1063         h_mesonClusterChi2->Fill(leadChi2);
1064         h_mesonClusterChi2->Fill(subChi2);
1065         }
1066         if (InRange(diphoton.M(), pi0BkgrLowRange) || InRange(diphoton.M(), pi0BkgrHighRange)) {
1067         /* vec = pi0Bkgr; */
1068         }
1069         if (InRange(diphoton.M(), etaBkgrLowRange) || InRange(diphoton.M(), etaBkgrHighRange)) {
1070         /* vec = etaBkgr; */
1071         }
1072
1073         /* if (vec) { */
1074         /* vec->push_back(d); */
1075         /* } */
1076     }
1077     } // end loop over pairs of clusters
1078     h_nDiphotons->Fill(nDiphotons);
1079     h_nRecoPi0s->Fill(nRecoPi0s);
1080     h_nRecoEtas->Fill(nRecoEtas);
1081     return;
1082 }
1083
1084 /* void neutralMesonTSSA::FillPhiHists(std::string which, int index) */
1085 /* { */
1086 /*     /1* std::vector<Diphoton>* vec = nullptr; *1/ */
1087 /*     PhiHists* hists = nullptr; */
1088 /*     PhiHists* hists_lowEta = nullptr; */
1089 /*     PhiHists* hists_highEta = nullptr; */
1090
1091 /*     if (which == "pi0") { */
1092 /*  /1* vec = pi0s; *1/ */
1093 /*  hists = pi0Hists; */
1094 /*  hists_lowEta = pi0Hists_lowEta; */
1095 /*  hists_highEta = pi0Hists_highEta; */
1096 /*     } */
1097 /*     if (which == "eta") { */
1098 /*  /1* vec = etas; *1/ */
1099 /*  hists = etaHists; */
1100 /*  hists_lowEta = etaHists_lowEta; */
1101 /*  hists_highEta = etaHists_highEta; */
1102 /*     } */
1103 /*     if (which == "pi0bkgr") { */
1104 /*  /1* vec = pi0Bkgr; *1/ */
1105 /*  hists = pi0BkgrHists; */
1106 /*  hists_lowEta = pi0BkgrHists_lowEta; */
1107 /*  hists_highEta = pi0BkgrHists_highEta; */
1108 /*     } */
1109 /*     if (which == "etabkgr") { */
1110 /*  /1* vec = etaBkgr; *1/ */
1111 /*  hists = etaBkgrHists; */
1112 /*  hists_lowEta = etaBkgrHists_lowEta; */
1113 /*  hists_highEta = etaBkgrHists_highEta; */
1114 /*     } */
1115
1116 /*     /1* if (!vec || !hists || !hists_lowEta) { *1/ */
1117 /*     if (!hists || !hists_lowEta || !hists_highEta) { */
1118 /*  std::cout << PHWHERE << ":: Invalid arguments!" << std::endl; */
1119 /*  return; */
1120 /*     } */
1121
1122 /*     Diphoton d = vec->at(index); */
1123 /*     // phi and xF need beam-dependent considerations */
1124 /*     float phi = d.phi; */
1125 /*     // rotate phi away from global coordinates, into PHENIX coordinate system */
1126 /*     float phiblue = 999.9; */
1127 /*     float phiyellow = 999.9; */
1128 /*     phiblue = phi - PI/2.0; */
1129 /*     if (phiblue < -PI) phiblue += 2.0*PI; */
1130 /*     phiyellow = phi + PI/2.0; */
1131 /*     if (phiyellow > PI) phiyellow -= 2.0*PI; */
1132 /*     // xF is defined such that forward is to the north. For yellow beam, */
1133 /*     // invert this such that xF > 0 <--> forward direction */
1134 /*     float xFblue = d.xF; */
1135 /*     float xFyellow = -1.0*d.xF; */
1136
1137 /*     hists->phi_pT->FillHists(d.pT, phi); */
1138 /*     hists->phi_xF->FillHists(d.xF, phi); */
1139 /*     hists->phi_eta->FillHists(d.eta, phi); */
1140 /*     hists->phi_vtxz->FillHists(d.vtxz, phi); */
1141 /*     if (bspin == 1) { */
1142 /*  hists->phi_pT_blue_up->FillHists(d.pT, phiblue); */
1143 /*  hists->phi_xF_blue_up->FillHists(xFblue, phiblue); */
1144 /*  hists->phi_eta_blue_up->FillHists(d.eta, phiblue); */
1145 /*  hists->phi_vtxz_blue_up->FillHists(d.vtxz, phiblue); */
1146 /*  if (abs(d.eta) < 0.35) */
1147 /*  { */
1148 /*      hists_lowEta->phi_pT_blue_up->FillHists(d.pT, phiblue); */
1149 /*      hists_lowEta->phi_xF_blue_up->FillHists(xFblue, phiblue); */
1150 /*      hists_lowEta->phi_eta_blue_up->FillHists(d.eta, phiblue); */
1151 /*      hists_lowEta->phi_vtxz_blue_up->FillHists(d.vtxz, phiblue); */
1152 /*  } */
1153 /*  if (abs(d.eta) > 0.35) */
1154 /*  { */
1155 /*      hists_highEta->phi_pT_blue_up->FillHists(d.pT, phiblue); */
1156 /*      hists_highEta->phi_xF_blue_up->FillHists(xFblue, phiblue); */
1157 /*      hists_highEta->phi_eta_blue_up->FillHists(d.eta, phiblue); */
1158 /*      hists_highEta->phi_vtxz_blue_up->FillHists(d.vtxz, phiblue); */
1159 /*  } */
1160 /*     } */
1161 /*     if (bspin == -1) { */
1162 /*  hists->phi_pT_blue_down->FillHists(d.pT, phiblue); */
1163 /*  hists->phi_xF_blue_down->FillHists(xFblue, phiblue); */
1164 /*  hists->phi_eta_blue_down->FillHists(d.eta, phiblue); */
1165 /*  hists->phi_vtxz_blue_down->FillHists(d.vtxz, phiblue); */
1166 /*  if (abs(d.eta) < 0.35) */
1167 /*  { */
1168 /*      hists_lowEta->phi_pT_blue_down->FillHists(d.pT, phiblue); */
1169 /*      hists_lowEta->phi_xF_blue_down->FillHists(xFblue, phiblue); */
1170 /*      hists_lowEta->phi_eta_blue_down->FillHists(d.eta, phiblue); */
1171 /*      hists_lowEta->phi_vtxz_blue_down->FillHists(d.vtxz, phiblue); */
1172 /*  } */
1173 /*  if (abs(d.eta) > 0.35) */
1174 /*  { */
1175 /*      hists_highEta->phi_pT_blue_down->FillHists(d.pT, phiblue); */
1176 /*      hists_highEta->phi_xF_blue_down->FillHists(xFblue, phiblue); */
1177 /*      hists_highEta->phi_eta_blue_down->FillHists(d.eta, phiblue); */
1178 /*      hists_highEta->phi_vtxz_blue_down->FillHists(d.vtxz, phiblue); */
1179 /*  } */
1180 /*     } */
1181 /*     if (yspin == 1) { */
1182 /*  hists->phi_pT_yellow_up->FillHists(d.pT, phiyellow); */
1183 /*  hists->phi_xF_yellow_up->FillHists(xFyellow, phiyellow); */
1184 /*  hists->phi_eta_yellow_up->FillHists(d.eta, phiyellow); */
1185 /*  hists->phi_vtxz_yellow_up->FillHists(d.vtxz, phiyellow); */
1186 /*  if (abs(d.eta) < 0.35) */
1187 /*  { */
1188 /*      hists_lowEta->phi_pT_yellow_up->FillHists(d.pT, phiyellow); */
1189 /*      hists_lowEta->phi_xF_yellow_up->FillHists(xFyellow, phiyellow); */
1190 /*      hists_lowEta->phi_eta_yellow_up->FillHists(d.eta, phiyellow); */
1191 /*      hists_lowEta->phi_vtxz_yellow_up->FillHists(d.vtxz, phiyellow); */
1192 /*  } */
1193 /*  if (abs(d.eta) > 0.35) */
1194 /*  { */
1195 /*      hists_highEta->phi_pT_yellow_up->FillHists(d.pT, phiyellow); */
1196 /*      hists_highEta->phi_xF_yellow_up->FillHists(xFyellow, phiyellow); */
1197 /*      hists_highEta->phi_eta_yellow_up->FillHists(d.eta, phiyellow); */
1198 /*      hists_highEta->phi_vtxz_yellow_up->FillHists(d.vtxz, phiyellow); */
1199 /*  } */
1200 /*     } */
1201 /*     if (yspin == -1) { */
1202 /*  hists->phi_pT_yellow_down->FillHists(d.pT, phiyellow); */
1203 /*  hists->phi_xF_yellow_down->FillHists(xFyellow, phiyellow); */
1204 /*  hists->phi_eta_yellow_down->FillHists(d.eta, phiyellow); */
1205 /*  hists->phi_vtxz_yellow_down->FillHists(d.vtxz, phiyellow); */
1206 /*  if (abs(d.eta) < 0.35) */
1207 /*  { */
1208 /*      hists_lowEta->phi_pT_yellow_down->FillHists(d.pT, phiyellow); */
1209 /*      hists_lowEta->phi_xF_yellow_down->FillHists(xFyellow, phiyellow); */
1210 /*      hists_lowEta->phi_eta_yellow_down->FillHists(d.eta, phiyellow); */
1211 /*      hists_lowEta->phi_vtxz_yellow_down->FillHists(d.vtxz, phiyellow); */
1212 /*  } */
1213 /*  if (abs(d.eta) > 0.35) */
1214 /*  { */
1215 /*      hists_highEta->phi_pT_yellow_down->FillHists(d.pT, phiyellow); */
1216 /*      hists_highEta->phi_xF_yellow_down->FillHists(xFyellow, phiyellow); */
1217 /*      hists_highEta->phi_eta_yellow_down->FillHists(d.eta, phiyellow); */
1218 /*      hists_highEta->phi_vtxz_yellow_down->FillHists(d.vtxz, phiyellow); */
1219 /*  } */
1220 /*     } */
1221 /* } */
1222
1223 /* void neutralMesonTSSA::FillAllPhiHists() */
1224 /* { */
1225 /*     for (unsigned int i=0; i<pi0s->size(); i++) { */
1226 /*  FillPhiHists("pi0", i); */
1227 /*     } */
1228 /*     for (unsigned int i=0; i<etas->size(); i++) { */
1229 /*  FillPhiHists("eta", i); */
1230 /*     } */
1231 /*     for (unsigned int i=0; i<pi0Bkgr->size(); i++) { */
1232 /*  FillPhiHists("pi0bkgr", i); */
1233 /*     } */
1234 /*     for (unsigned int i=0; i<etaBkgr->size(); i++) { */
1235 /*  FillPhiHists("etabkgr", i); */
1236 /*     } */
1237 /* } */
1238
1239 void neutralMesonTSSA::ClearVectors()
1240 {
1241     goodclusters_E->clear();
1242     goodclusters_Ecore->clear();
1243     goodclusters_Eta->clear();
1244     goodclusters_Phi->clear();
1245     goodclusters_pT->clear();
1246     goodclusters_xF->clear();
1247     goodclusters_Chi2->clear();
1248
1249     diphoton_E->clear();
1250     diphoton_M->clear();
1251     diphoton_Eta->clear();
1252     diphoton_Phi->clear();
1253     diphoton_pT->clear();
1254     diphoton_xF->clear();
1255     diphoton_clus1index->clear();
1256     diphoton_clus2index->clear();
1257     diphoton_deltaR->clear();
1258     diphoton_asym->clear();
1259
1260     /* pi0s->clear(); */
1261     /* etas->clear(); */
1262     /* pi0Bkgr->clear(); */
1263     /* etaBkgr->clear(); */
1264 }
1265
1266 void neutralMesonTSSA::DeleteStuff()
1267 {
1268     /* delete pi0Hists; delete etaHists; delete pi0BkgrHists; delete etaBkgrHists; */
1269     delete goodclusters_E; delete goodclusters_Ecore; delete goodclusters_Eta; delete goodclusters_Phi; delete goodclusters_pT; delete goodclusters_xF; delete goodclusters_Chi2;
1270     delete diphoton_E; delete diphoton_M; delete diphoton_Eta; delete diphoton_Phi; delete diphoton_pT; delete diphoton_xF; delete diphoton_clus1index; delete diphoton_clus2index; delete diphoton_deltaR; delete diphoton_asym;
1271     /* delete pi0s; delete etas; delete pi0Bkgr; delete etaBkgr; */
1272     delete outfile_hists; delete outfile_trees;
1273 }
1274
1275 void neutralMesonTSSA::GetTrigger()
1276 {
1277     if (gl1Packet)
1278     {
1279     /* uint64_t trig = gl1Packet->getTriggerVector(); */
1280     uint64_t live = gl1Packet->getLiveVector();
1281     uint64_t scaled = gl1Packet->getScaledVector();
1282
1283     unsigned int mbdtrignum = 10;
1284     if (((live >> mbdtrignum ) & 0x1U ) == 0x1U)
1285     {
1286         mbdtrigger_live = true;
1287         /* n_events_mbdtrigger++; */
1288     }
1289     if (((scaled >> mbdtrignum ) & 0x1U ) == 0x1U)
1290     {
1291         mbdtrigger_scaled = true;
1292         n_events_mbdtrigger++;
1293     }
1294
1295     unsigned int photontrignum = 25; // 3GeV photon + MBD coinc. trigger
1296     if (((live >> photontrignum ) & 0x1U ) == 0x1U)
1297     {
1298         photontrigger_live = true;
1299         /* n_events_photontrigger++; */
1300     }
1301     if (((scaled >> photontrignum ) & 0x1U ) == 0x1U)
1302     {
1303         photontrigger_scaled = true;
1304         n_events_photontrigger++;
1305     }
1306     }
1307 }
1308
1309 void neutralMesonTSSA::PrintTrigger()
1310 {
1311     if (gl1Packet)
1312     {
1313     uint64_t trig = gl1Packet->getTriggerVector();
1314     uint64_t live = gl1Packet->getLiveVector();
1315     uint64_t scaled = gl1Packet->getScaledVector();
1316
1317     std::bitset<64> trigbits(trig);
1318     std::bitset<64> livebits(live);
1319     std::bitset<64> scaledbits(scaled);
1320
1321     /* std::cout << "trig = " << trig << "\ntrigbits = " << trigbits << std::endl; */
1322     for (unsigned int i=0; i<64; i++)
1323     {
1324         if (((trig >> i ) & 0x1U ) == 0x1U)
1325         {
1326         /* std::cout << "Trigger: bit " << i << " = true" << std::endl; */
1327         if (i == 8) mbdStrigger = true;
1328         if (i == 9) mbdNtrigger = true;
1329         if (i == 10) {
1330             n_events_mbdtrigger++;
1331             mbdtrigger_live = true;
1332             if (!(mbdStrigger && mbdNtrigger)) mbdcoinc_withoutNandS++;
1333         }
1334         if (i == 12) n_events_mbdtrigger_vtx1++;
1335         if (i == 13) n_events_mbdtrigger_vtx2++;
1336         if (i == 14) n_events_mbdtrigger_vtx3++;
1337         }
1338     }
1339     /*
1340     std::cout << "live = " << live << "\nlivebits = " << livebits << std::endl;
1341     for (unsigned int i=0; i<64; i++)
1342     {
1343         if (((live >> i ) & 0x1U ) == 0x1U)
1344         {
1345         std::cout << "Live: bit " << i << " = true" << std::endl;
1346         }
1347     }
1348     */
1349     /*
1350     std::cout << "scaled = " << scaled << "\nscaledbits = " << scaledbits << std::endl;
1351     for (unsigned int i=0; i<64; i++)
1352     {
1353         if (((scaled >> i ) & 0x1U ) == 0x1U)
1354         {
1355         std::cout << "Scaled: bit " << i << " = true" << std::endl;
1356         }
1357     }
1358     */
1359     }
1360 }