TSSAhistmaker/src/TSSAhistmaker.cc

0001 #include "TSSAhistmaker.h"
0002 #include "binnedhistset/BinnedHistSet.h"
0003
0004 #include <TFile.h>
0005 #include <TTree.h>
0006 #include <TChain.h>
0007 #include <TLorentzVector.h>
0008 #include <TVector3.h>
0009 #include <TH1.h>
0010 #include <TString.h>
0011
0012 #include <iostream>
0013
0014 TSSAhistmaker::TSSAhistmaker(const std::string treepath, const int jobnumber, const std::string histfname):
0015     treebasepath(treepath),
0016     jobnum(jobnumber),
0017     histfilename(histfname)
0018 {
0019 }
0020
0021 TSSAhistmaker::~TSSAhistmaker() {}
0022
0023 bool TSSAhistmaker::InRange(float mass, std::pair<float, float> range)
0024 {
0025     bool ret = false;
0026     if ( (mass > range.first) && (mass < range.second) ) ret = true;
0027     return ret;
0028 }
0029
0030 void TSSAhistmaker::GetTrees()
0031 {
0032     /* infile_trees = new TFile(treefilename.c_str(), "READ"); */
0033     /* tree_clusters = (TTree*)infile_trees->Get("Clusters"); */
0034     /* tree_diphotons = (TTree*)infile_trees->Get("Diphotons"); */
0035     /* if (!tree_clusters) { */
0036     /* std::cout << "Unable to read tree Clusters from " << treefilename << "; exiting!"; */
0037     /* exit(1); */
0038     /* } */
0039     /* if (!tree_diphotons) { */
0040     /* std::cout << "Unable to read tree Diphotons from " << treefilename << "; exiting!"; */
0041     /* exit(1); */
0042     /* } */
0043     tree_event = new TChain("Event");
0044     tree_clusters = new TChain("Clusters");
0045     tree_diphotons = new TChain("Diphotons");
0046     int startingtree = jobnum*nTreesInChain;
0047     for (int i=startingtree; i<(startingtree+nTreesInChain); i++) {
0048     std::string treepath = treebasepath + "/" + std::to_string(i) + "/";
0049     treefilename = treepath + std::string("NMtrees-") + std::to_string(i) + std::string(".root");
0050     std::cout << "Adding " << treefilename << " to TChains" << std::endl;;
0051     tree_event->Add(treefilename.c_str());
0052     tree_clusters->Add(treefilename.c_str());
0053     tree_diphotons->Add(treefilename.c_str());
0054     /* std::cout << "Adding " << treefilename << " to TChains... "; */
0055     /* int cluschain = tree_clusters->AddFile(treefilename.c_str(), TTree::kMaxEntries, "Clusters"); */
0056     /* int diphchain = tree_diphotons->AddFile(treefilename.c_str(), TTree::kMaxEntries, "Diphotons"); */
0057     /* std::cout << cluschain << ", " << diphchain << std::endl; */
0058     }
0059
0060     /* tree_event->Print(); */
0061     /* tree_clusters->Print(); */
0062     /* tree_diphotons->Print(); */
0063     tree_event->SetBranchAddress("crossingAngle", &crossingAngle);
0064     tree_event->SetBranchAddress("vtxz", &vtxz);
0065     tree_event->SetBranchAddress("minbiastrig_live", &minbiastrig_live);
0066     tree_event->SetBranchAddress("photontrig_live", &photontrig_live);
0067     tree_event->SetBranchAddress("minbiastrig_scaled", &minbiastrig_scaled);
0068     tree_event->SetBranchAddress("photontrig_scaled", &photontrig_scaled);
0069     tree_event->SetBranchAddress("bspin", &bspin);
0070     tree_event->SetBranchAddress("yspin", &yspin);
0071
0072     tree_clusters->SetBranchAddress("clusterE", &cluster_E);
0073     tree_clusters->SetBranchAddress("clusterEcore", &cluster_Ecore);
0074     tree_clusters->SetBranchAddress("clusterEta", &cluster_Eta);
0075     tree_clusters->SetBranchAddress("clusterPhi", &cluster_Phi);
0076     tree_clusters->SetBranchAddress("clusterpT", &cluster_pT);
0077     tree_clusters->SetBranchAddress("clusterxF", &cluster_xF);
0078     tree_clusters->SetBranchAddress("clusterChi2", &cluster_Chi2);
0079     tree_clusters->SetBranchAddress("nAllClusters", &nAllClusters);
0080
0081     tree_diphotons->SetBranchAddress("diphotonE", &diphoton_E);
0082     tree_diphotons->SetBranchAddress("diphotonM", &diphoton_M);
0083     tree_diphotons->SetBranchAddress("diphotonEta", &diphoton_Eta);
0084     tree_diphotons->SetBranchAddress("diphotonPhi", &diphoton_Phi);
0085     tree_diphotons->SetBranchAddress("diphotonpT", &diphoton_pT);
0086     tree_diphotons->SetBranchAddress("diphotonxF", &diphoton_xF);
0087     tree_diphotons->SetBranchAddress("diphotonClusterIndex1", &diphoton_clus1index);
0088     tree_diphotons->SetBranchAddress("diphotonClusterIndex2", &diphoton_clus2index);
0089     tree_diphotons->SetBranchAddress("diphotonDeltaR", &diphoton_deltaR);
0090     tree_diphotons->SetBranchAddress("diphotonAsym", &diphoton_asym);
0091 }
0092
0093 void TSSAhistmaker::MakePhiHists(std::string which)
0094 {
0095     PhiHists* hists = new PhiHists();
0096     std::string nameprefix, titlewhich, titlebeam;
0097     nameprefix = "h_" + which + "_";
0098     std::vector<double> pTbins;
0099
0100     if (which == "pi0") {
0101     titlewhich = "#pi^{0}";
0102     pi0Hists = hists;
0103     pTbins = pTbins_pi0;
0104     }
0105     else if (which == "eta") {
0106     titlewhich = "#eta";
0107     etaHists = hists;
0108     pTbins = pTbins_eta;
0109     }
0110     else if (which == "pi0bkgr") {
0111     titlewhich = "#pi^{0} Background";
0112     pi0BkgrHists = hists;
0113     pTbins = pTbins_pi0;
0114     }
0115     else if (which == "etabkgr") {
0116     titlewhich = "#eta Background";
0117     etaBkgrHists = hists;
0118     pTbins = pTbins_eta;
0119     }
0120     else if (which == "pi0_lowEta") {
0121     titlewhich = "#pi^{0}, |#eta| < 0.35";
0122     pi0Hists_lowEta = hists;
0123     pTbins = pTbins_pi0;
0124     }
0125     else if (which == "eta_lowEta") {
0126     titlewhich = "#eta, |#eta| < 0.35";
0127     etaHists_lowEta = hists;
0128     pTbins = pTbins_eta;
0129     }
0130     else if (which == "pi0bkgr_lowEta") {
0131     titlewhich = "#pi^{0} Background, |#eta| < 0.35";
0132     pi0BkgrHists_lowEta = hists;
0133     pTbins = pTbins_pi0;
0134     }
0135     else if (which == "etabkgr_lowEta") {
0136     titlewhich = "#eta Background, |#eta| < 0.35";
0137     etaBkgrHists_lowEta = hists;
0138     pTbins = pTbins_eta;
0139     }
0140     else if (which == "pi0_highEta") {
0141     titlewhich = "#pi^{0}, |#eta| > 0.35";
0142     pi0Hists_highEta = hists;
0143     pTbins = pTbins_pi0;
0144     }
0145     else if (which == "eta_highEta") {
0146     titlewhich = "#eta, |#eta| > 0.35";
0147     etaHists_highEta = hists;
0148     pTbins = pTbins_eta;
0149     }
0150     else if (which == "pi0bkgr_highEta") {
0151     titlewhich = "#pi^{0} Background, |#eta| > 0.35";
0152     pi0BkgrHists_highEta = hists;
0153     pTbins = pTbins_pi0;
0154     }
0155     else if (which == "etabkgr_highEta") {
0156     titlewhich = "#eta Background, |#eta| > 0.35";
0157     etaBkgrHists_highEta = hists;
0158     pTbins = pTbins_eta;
0159     }
0160     else {
0161     std::cout << "In MakePhiHists(): Invalid argument " << which << " !" << std::endl;
0162     return;
0163     }
0164
0165     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);
0166     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);
0167     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);
0168     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);
0169     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);
0170     hists->phi_pT_blue_up_forward = new BinnedHistSet(Form("%sphi_pT_blue_up_forward", nameprefix.c_str()), Form("%s Blue Beam Forward-Going Spin-Up #phi^{B} Distribution;#phi^{B} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "p_{T} (GeV)", pTbins);
0171     hists->phi_pT_blue_down_forward = new BinnedHistSet(Form("%sphi_pT_blue_down_forward", nameprefix.c_str()), Form("%s Blue Beam Forward-Going Spin-Down #phi^{B} Distribution;#phi^{B} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "p_{T} (GeV)", pTbins);
0172     hists->phi_pT_yellow_up_forward = new BinnedHistSet(Form("%sphi_pT_yellow_up_forward", nameprefix.c_str()), Form("%s Yellow Beam Forward-Going Spin-Up #phi^{Y} Distribution;#phi^{Y} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "p_{T} (GeV)", pTbins);
0173     hists->phi_pT_yellow_down_forward = new BinnedHistSet(Form("%sphi_pT_yellow_down_forward", nameprefix.c_str()), Form("%s Yellow Beam Forward-Going Spin-Down #phi^{Y} Distribution;#phi^{Y} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "p_{T} (GeV)", pTbins);
0174     hists->phi_pT_blue_up_backward = new BinnedHistSet(Form("%sphi_pT_blue_up_backward", nameprefix.c_str()), Form("%s Blue Beam Backward-Going Spin-Up #phi^{B} Distribution;#phi^{B} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "p_{T} (GeV)", pTbins);
0175     hists->phi_pT_blue_down_backward = new BinnedHistSet(Form("%sphi_pT_blue_down_backward", nameprefix.c_str()), Form("%s Blue Beam Backward-Going Spin-Down #phi^{B} Distribution;#phi^{B} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "p_{T} (GeV)", pTbins);
0176     hists->phi_pT_yellow_up_backward = new BinnedHistSet(Form("%sphi_pT_yellow_up_backward", nameprefix.c_str()), Form("%s Yellow Beam Backward-Going Spin-Up #phi^{Y} Distribution;#phi^{Y} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "p_{T} (GeV)", pTbins);
0177     hists->phi_pT_yellow_down_backward = new BinnedHistSet(Form("%sphi_pT_yellow_down_backward", nameprefix.c_str()), Form("%s Yellow Beam Backward-Going Spin-Down #phi^{Y} Distribution;#phi^{Y} (rad);Counts", titlewhich.c_str()), nHistBins_phi, -1.0*PI, PI, "p_{T} (GeV)", pTbins);
0178
0179     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);
0180     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);
0181     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);
0182     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);
0183     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);
0184
0185     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);
0186     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);
0187     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);
0188     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);
0189     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);
0190
0191     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);
0192     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);
0193     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);
0194     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);
0195     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);
0196
0197     hists->phi_pT->MakeHists();
0198     hists->phi_pT_blue_up->MakeHists();
0199     hists->phi_pT_blue_down->MakeHists();
0200     hists->phi_pT_yellow_up->MakeHists();
0201     hists->phi_pT_yellow_down->MakeHists();
0202     hists->phi_pT_blue_up_forward->MakeHists();
0203     hists->phi_pT_blue_down_forward->MakeHists();
0204     hists->phi_pT_yellow_up_forward->MakeHists();
0205     hists->phi_pT_yellow_down_forward->MakeHists();
0206     hists->phi_pT_blue_up_backward->MakeHists();
0207     hists->phi_pT_blue_down_backward->MakeHists();
0208     hists->phi_pT_yellow_up_backward->MakeHists();
0209     hists->phi_pT_yellow_down_backward->MakeHists();
0210
0211     hists->phi_xF->MakeHists();
0212     hists->phi_xF_blue_up->MakeHists();
0213     hists->phi_xF_blue_down->MakeHists();
0214     hists->phi_xF_yellow_up->MakeHists();
0215     hists->phi_xF_yellow_down->MakeHists();
0216
0217     hists->phi_eta->MakeHists();
0218     hists->phi_eta_blue_up->MakeHists();
0219     hists->phi_eta_blue_down->MakeHists();
0220     hists->phi_eta_yellow_up->MakeHists();
0221     hists->phi_eta_yellow_down->MakeHists();
0222
0223     hists->phi_vtxz->MakeHists();
0224     hists->phi_vtxz_blue_up->MakeHists();
0225     hists->phi_vtxz_blue_down->MakeHists();
0226     hists->phi_vtxz_yellow_up->MakeHists();
0227     hists->phi_vtxz_yellow_down->MakeHists();
0228 }
0229
0230 void TSSAhistmaker::MakeAllHists()
0231 {
0232     // create file
0233     outfile_hists = new TFile(histfilename.c_str(), "RECREATE");
0234     outfile_hists->cd();
0235
0236     // event-level
0237     int nbins_vtxz = 100;
0238     double vtxz_upper = 200.0;
0239     h_vtxz = new TH1F("h_vtxz", "Vertex z Distribution;z_{vtx} (cm);Counts", nbins_vtxz, -vtxz_upper, vtxz_upper);
0240     h_Nevents = new TH1F("h_Nevents", "Number of Events (Minbias trig, GlobalVertex, Positive E_{EMCal}, N_{clusters} > 1);;Counts", 1, -0.5, 0.5);
0241     h_Npi0s = new TH1F("h_Npi0s", "Number of Diphotons in #pi^{0} Mass Range;;Counts", 1, -0.5, 0.5);
0242     h_Netas = new TH1F("h_Netas", "Number of Diphotons in #eta Mass Range;;Counts", 1, -0.5, 0.5);
0243
0244     // clusters
0245     /* h_nAllClusters = new TH1F("h_nAllClusters", "Total Number of Clusters per Event;# Clusters;Counts", 3500, 0.0, 3500.0); */
0246     /* h_nGoodClusters = new TH1F("h_nGoodClusters", "Number of \"Good\" Clusters per Event;# Clusters;Counts", 13, -0.5, 12.5); */
0247     h_nClustersMinbiasTrig = new TH1F("h_nClustersMinbiasTrig", "Number of \"Good\" Clusters per Event, Minbias Trigger;# Clusters;Counts", 13, -0.5, 12.5);
0248     h_nClustersPhotonTrig = new TH1F("h_nClustersPhotonTrig", "Number of \"Good\" Clusters per Event, Photon Trigger;# Clusters;Counts", 13, -0.5, 12.5);
0249
0250     // diphotons
0251     int nbins_mass = 100;
0252     int nbins_asym = 100;
0253     int nbins_dR = 100;
0254     double dR_upper_small = 0.5;
0255     h_DiphotonMassAsym = new TH2F("h_DiphotonMassAsym", "Diphoton Pair Asymmetry vs Mass;Mass (GeV);#alpha", nbins_mass, 0.0, 1.0, nbins_asym, 0.0, 1.0);
0256     h_DiphotonMassdeltaR = new TH2F("h_DiphotonMassdeltaR", "Diphoton #Delta R vs Mass; Mass (GeV);#Delta R", nbins_mass, 0.0, 1.0, nbins_dR, 0.0, PI);
0257     h_DiphotondeltaRAsym = new TH2F("h_DiphotondeltaRAsym", "Diphoton Pair Asymmetry vs #Delta R;#Delta R;#alpha", nbins_dR, 0.0, PI, nbins_dR, 0.0, 1.0);
0258     h_DiphotondeltaRAsym_pi0 = new TH2F("h_DiphotondeltaRAsym_pi0", "Diphoton Pair Asymmetry vs #Delta R, Mass in #pi^{0} Range;#Delta R;#alpha", nbins_dR, 0.0, PI, nbins_dR, 0.0, 1.0);
0259     h_DiphotondeltaRAsym_pi0_smalldR = new TH2F("h_DiphotondeltaRAsym_pi0_smalldR", "Diphoton Pair Asymmetry vs #Delta R, Mass in #pi^{0} Range;#Delta R;#alpha", nbins_dR, 0.0, dR_upper_small, nbins_dR, 0.0, 1.0);
0260     h_DiphotondeltaRAsym_eta = new TH2F("h_DiphotondeltaRAsym_eta", "Diphoton Pair Asymmetry vs #Delta R, Mass in #eta Range;#Delta R;#alpha", nbins_dR, 0.0, PI, nbins_dR, 0.0, 1.0);
0261     h_DiphotondeltaRAsym_etabkgr = new TH2F("h_DiphotondeltaRAsym_etabkgr", "Diphoton Pair Asymmetry vs #Delta R, Mass in #eta Background Range;#Delta R;#alpha", nbins_dR, 0.0, PI, nbins_dR, 0.0, 1.0);
0262     h_DiphotonMassAsym_highpT = new TH2F("h_DiphotonMassAsym_highpT", "Diphoton Pair Asymmetry vs Mass, p_{T} > 7 GeV;Mass (GeV);#alpha", nbins_mass, 0.0, 1.0, nbins_asym, 0.0, 1.0);
0263     h_DiphotonMassdeltaR_highpT = new TH2F("h_DiphotonMassdeltaR_highpT", "Diphoton #Delta R vs Mass, p_{T} > 7 GeV;Mass (GeV);#Delta R", nbins_mass, 0.0, 1.0, nbins_dR, 0.0, dR_upper_small);
0264     h_DiphotondeltaRAsym_highpT = new TH2F("h_DiphotondeltaRAsym_highpT", "Diphoton Pair Asymmetry vs #Delta R, p_{T} > 7 GeV;#Delta R;#alpha", nbins_dR, 0.0, PI, nbins_dR, 0.0, 1.0);
0265     h_DiphotondeltaRAsym_highpT_smalldR = new TH2F("h_DiphotondeltaRAsym_highpT_smalldR", "Diphoton Pair Asymmetry vs #Delta R, p_{T} > 7 GeV;#Delta R;#alpha", nbins_dR, 0.0, dR_upper_small, nbins_dR, 0.0, 1.0);
0266     h_DiphotondeltaRAsym_eta_highpT = new TH2F("h_DiphotondeltaRAsym_eta_highpT", "Diphoton Pair Asymmetry vs #Delta R, Mass in #eta Range, p_{T} > 7 GeV;#Delta R;#alpha", nbins_dR, 0.0, dR_upper_small, nbins_dR, 0.0, 1.0);
0267     h_DiphotondeltaRAsym_etabkgr_highpT = new TH2F("h_DiphotondeltaRAsym_etabkgr_highpT", "Diphoton Pair Asymmetry vs #Delta R, Mass in #eta Background Range, p_{T} > 7 GeV;#Delta R;#alpha", nbins_dR, 0.0, dR_upper_small, nbins_dR, 0.0, 1.0);
0268     // Armenteros-Podolanski plots
0269     h_armen_all = new TH2F("h_armen_all", "Armenteros-Podolanski Plot, All Diphotons;#alpha = #frac{p_{L1} - p_{L2}}{p_{L1} + p_{L2}};p_{T_{#gamma#gamma}}", 100, -1.0, 1.0, 100, 0.0, 5.0);
0270     h_armen_pi0 = new TH2F("h_armen_pi0", "Armenteros-Podolanski Plot, #pi^{0} Mass Range;#alpha = #frac{p_{L1} - p_{L2}}{p_{L1} + p_{L2}};p_{T_{#gamma#gamma}}", 100, -1.0, 1.0, 100, 0.0, 0.5);
0271     h_armen_pi0bkgr = new TH2F("h_armen_pi0bkgr", "Armenteros-Podolanski Plot, #pi^{0} Background Mass Range;#alpha = #frac{p_{L1} - p_{L2}}{p_{L1} + p_{L2}};p_{T_{#gamma#gamma}}", 100, -1.0, 1.0, 100, 0.0, 0.5);
0272     h_armen_eta = new TH2F("h_armen_eta", "Armenteros-Podolanski Plot, #eta Mass Range;#alpha = #frac{p_{L1} - p_{L2}}{p_{L1} + p_{L2}};p_{T_{#gamma#gamma}}", 100, -1.0, 1.0, 100, 0.0, 0.5);
0273     h_armen_etabkgr = new TH2F("h_armen_etabkgr", "Armenteros-Podolanski Plot, #eta Background Mass Range;#alpha = #frac{p_{L1} - p_{L2}}{p_{L1} + p_{L2}};p_{T_{#gamma#gamma}}", 100, -1.0, 1.0, 100, 0.0, 0.5);
0274     h_armen_eta_highpT = new TH2F("h_armen_eta_highpT", "Armenteros-Podolanski Plot, #eta Mass Range, p_{T} > 7 GeV;#alpha = #frac{p_{L1} - p_{L2}}{p_{L1} + p_{L2}};p_{T_{#gamma#gamma}}", 100, -1.0, 1.0, 100, 0.0, 0.5);
0275     h_armen_etabkgr_highpT = new TH2F("h_armen_etabkgr_highpT", "Armenteros-Podolanski Plot, #eta Background Mass Range, p_{T} > 7 GeV;#alpha = #frac{p_{L1} - p_{L2}}{p_{L1} + p_{L2}};p_{T_{#gamma#gamma}}", 100, -1.0, 1.0, 100, 0.0, 0.5);
0276     h_armen_p_L = new TH1F("h_armen_p_L", "Cluster p_{L} w.r.t. Diphoton;p_{L, #gamma#gamma} (GeV);Counts", 200, 0.0, 10.0);
0277     h_armen_p_T = new TH1F("h_armen_p_T", "Cluster p_{T} w.r.t. Diphoton;p_{T, #gamma#gamma} (GeV);Counts", 200, 0.0, 10.0);
0278     h_armen_alpha_alphaE = new TH2F("h_armen_alpha_alphaE", "Energy Asymmetry vs Armenteros #alpha;(p_{L1} - p_{L2})/(p_{L1} + p_{L2});(E1 - E2)/(E1 + E2)", 100, -1.0, 1.0, 100, -1.0, 1.0);
0279     h_armen_alpha_deltaR = new TH2F("h_armen_alpha_deltaR", "Cluster #Delta R vs Armenteros #alpha;(p_{L1} - p_{L2})/(p_{L1} + p_{L2});#Delta R (rad)", 100, -1.0, 1.0, 100, 0.0, PI);
0280     h_armen_alpha_deltaR_pi0 = new TH2F("h_armen_alpha_deltaR_pi0", "Cluster #Delta R vs Armenteros #alpha, #pi^{0} Mass Range;(p_{L1} - p_{L2})/(p_{L1} + p_{L2});#Delta R (rad)", 100, -1.0, 1.0, 100, 0.0, PI);
0281     h_armen_alpha_deltaR_eta = new TH2F("h_armen_alpha_deltaR_eta", "Cluster #Delta R vs Armenteros #alpha, #eta Mass Range;(p_{L1} - p_{L2})/(p_{L1} + p_{L2});#Delta R (rad)", 100, -1.0, 1.0, 100, 0.0, PI);
0282     h_armen_p_T_deltaR = new TH2F("h_armen_p_T_deltaR", "Cluster #Delta R vs p_{T} w.r.t. Diphoton;p_{T, #gamma#gamma};#Delta R (rad)", 200, 0.0, 10.0, 100, 0.0, PI);
0283     h_armen_pT_p_L = new TH2F("h_armen_pT_p_L", "Diphoton p_{T} vs Cluster p_{L} w.r.t. Diphoton;p_{T} (GeV);p_{L, #gamma#gamma}", 200, 0.0, 20.0, 100, 0.0, 10.0);
0284     h_armen_pT_p_T = new TH2F("h_armen_pT_p_T", "Diphoton p_{T} vs Cluster p_{T} w.r.t. Diphoton;p_{T} (GeV);p_{T, #gamma#gamma}", 200, 0.0, 20.0, 100, 0.0, 10.0);
0285
0286     // diphoton mass
0287     float pT_upper = 20.0;
0288     h_diphotonMass = new TH1F("h_diphotonMass", "Diphoton Mass Distribution (Cluster p_{T} > 1.5GeV);Diphoton Mass (GeV);Counts", nbins_mass, 0.0, 1.0);
0289     std::vector<double> pTbins_mass;
0290     /* std::vector<double> xFbins; */
0291     int nbins_bhs = 20;
0292     for (int i=0; i<nbins_bhs; i++) {
0293     pTbins_mass.push_back(i*(pT_upper/nbins_bhs));
0294     /* xFbins.push_back(2*xF_upper*i/nbins_bhs - xF_upper); */
0295     }
0296     pTbins_mass.push_back(pT_upper);
0297     /* xFbins.push_back(xF_upper); */
0298     /* bhs_diphotonMass_pT = new BinnedHistSet("h_diphotonMass_pT", "Diphoton Mass;Mass (GeV);Counts", nbins_mass, 0.0, 1.0, "p_{T} (GeV)", pTbins_mass); */
0299     bhs_diphotonMass_pT_pi0binning = new BinnedHistSet("h_diphotonMass_pT_pi0binning", "Diphoton Mass;Mass (GeV);Counts", nbins_mass, 0.0, 1.0, "p_{T} (GeV)", pTbins_pi0);
0300     bhs_diphotonMass_pT_etabinning = new BinnedHistSet("h_diphotonMass_pT_etabinning", "Diphoton Mass;Mass (GeV);Counts", nbins_mass, 0.0, 1.0, "p_{T} (GeV)", pTbins_eta);
0301     bhs_pi0_pT_pT = new BinnedHistSet("h_pi0_pT_pT", "#pi^{0} p_{T};p_{T} (GeV);Counts", 100, 1.0, 10.0, "p_{T} (GeV)", pTbins_pi0);
0302     bhs_eta_pT_pT = new BinnedHistSet("h_eta_pT_pT", "#eta p_{T};p_{T} (GeV);Counts", 100, 2.0, 20.0, "p_{T} (GeV)", pTbins_eta);
0303     bhs_diphotonxF_xF = new BinnedHistSet("h_diphotonxF_xF", "Diphoton x_{F};x_{F};Counts", 100, -0.15, 0.15, "x_{F}", xFbins);
0304     bhs_diphotoneta_eta = new BinnedHistSet("h_diphotoneta_eta", "Diphoton #eta;#eta;Counts", 100, -3.0, 3.0, "#eta", etabins);
0305     bhs_diphotonvtxz_vtxz = new BinnedHistSet("h_diphotonvtxz_vtxz", "Diphoton vtxz;vtxz (cm);Counts", 100, -200.0, 200.0, "vtxz (cm)", vtxzbins);
0306     bhs_diphotonMass_pT_pi0binning->MakeHists();
0307     bhs_diphotonMass_pT_etabinning->MakeHists();
0308     bhs_pi0_pT_pT->MakeHists();
0309     bhs_eta_pT_pT->MakeHists();
0310     bhs_diphotonxF_xF->MakeHists();
0311     bhs_diphotoneta_eta->MakeHists();
0312     bhs_diphotonvtxz_vtxz->MakeHists();
0313
0314     // phi hists
0315     MakePhiHists("pi0");
0316     MakePhiHists("eta");
0317     MakePhiHists("pi0bkgr");
0318     MakePhiHists("etabkgr");
0319     MakePhiHists("pi0_lowEta");
0320     MakePhiHists("eta_lowEta");
0321     MakePhiHists("pi0bkgr_lowEta");
0322     MakePhiHists("etabkgr_lowEta");
0323     MakePhiHists("pi0_highEta");
0324     MakePhiHists("eta_highEta");
0325     MakePhiHists("pi0bkgr_highEta");
0326     MakePhiHists("etabkgr_highEta");
0327 }
0328
0329 void TSSAhistmaker::FillPhiHists(std::string which, int index)
0330 {
0331     PhiHists* hists = nullptr;
0332     PhiHists* hists_lowEta = nullptr;
0333     PhiHists* hists_highEta = nullptr;
0334
0335     if (which == "pi0") {
0336     hists = pi0Hists;
0337     hists_lowEta = pi0Hists_lowEta;
0338     hists_highEta = pi0Hists_highEta;
0339     }
0340     if (which == "eta") {
0341     hists = etaHists;
0342     hists_lowEta = etaHists_lowEta;
0343     hists_highEta = etaHists_highEta;
0344     }
0345     if (which == "pi0bkgr") {
0346     hists = pi0BkgrHists;
0347     hists_lowEta = pi0BkgrHists_lowEta;
0348     hists_highEta = pi0BkgrHists_highEta;
0349     }
0350     if (which == "etabkgr") {
0351     hists = etaBkgrHists;
0352     hists_lowEta = etaBkgrHists_lowEta;
0353     hists_highEta = etaBkgrHists_highEta;
0354     }
0355
0356     if (!hists || !hists_lowEta || !hists_highEta) {
0357     std::cout << "In FillPhiHists: Invalid argument" << which << " !" << std::endl;
0358     return;
0359     }
0360
0361     float eta = diphoton_Eta->at(index);
0362     float phi = diphoton_Phi->at(index);
0363     float pT = diphoton_pT->at(index);
0364     float xF = diphoton_xF->at(index);
0365     // phi and xF need beam-dependent considerations
0366     // rotate phi away from global coordinates, into PHENIX coordinate system
0367     float phiblue = 999.9;
0368     float phiyellow = 999.9;
0369     phiblue = phi - PI/2.0;
0370     if (phiblue < -PI) phiblue += 2.0*PI;
0371     phiyellow = phi + PI/2.0;
0372     if (phiyellow > PI) phiyellow -= 2.0*PI;
0373     /* // xF is defined such that forward is to the north. For yellow beam, */
0374     /* // invert this such that xF > 0 <--> forward direction */
0375     float xFblue = xF;
0376     float xFyellow = -1.0*xF;
0377     /* float xFyellow = xF; */
0378
0379     hists->phi_pT->FillHists(pT, phi);
0380     hists->phi_xF->FillHists(xF, phi);
0381     hists->phi_eta->FillHists(eta, phi);
0382     hists->phi_vtxz->FillHists(vtxz, phi);
0383     if (bspin == 1) {
0384     hists->phi_pT_blue_up->FillHists(pT, phiblue);
0385     if (eta > 0) hists->phi_pT_blue_up_forward->FillHists(pT, phiblue);
0386     if (eta < 0) hists->phi_pT_blue_up_backward->FillHists(pT, phiblue);
0387     hists->phi_xF_blue_up->FillHists(xFblue, phiblue);
0388     hists->phi_eta_blue_up->FillHists(eta, phiblue);
0389     hists->phi_vtxz_blue_up->FillHists(vtxz, phiblue);
0390     if (abs(eta) < 0.35)
0391     {
0392         hists_lowEta->phi_pT_blue_up->FillHists(pT, phiblue);
0393         if (eta > 0) hists_lowEta->phi_pT_blue_up_forward->FillHists(pT, phiblue);
0394         if (eta < 0) hists_lowEta->phi_pT_blue_up_backward->FillHists(pT, phiblue);
0395         hists_lowEta->phi_xF_blue_up->FillHists(xFblue, phiblue);
0396         hists_lowEta->phi_eta_blue_up->FillHists(eta, phiblue);
0397         hists_lowEta->phi_vtxz_blue_up->FillHists(vtxz, phiblue);
0398     }
0399     if (abs(eta) > 0.35)
0400     {
0401         hists_highEta->phi_pT_blue_up->FillHists(pT, phiblue);
0402         if (eta > 0) hists_highEta->phi_pT_blue_up_forward->FillHists(pT, phiblue);
0403         if (eta < 0) hists_highEta->phi_pT_blue_up_backward->FillHists(pT, phiblue);
0404         hists_highEta->phi_xF_blue_up->FillHists(xFblue, phiblue);
0405         hists_highEta->phi_eta_blue_up->FillHists(eta, phiblue);
0406         hists_highEta->phi_vtxz_blue_up->FillHists(vtxz, phiblue);
0407     }
0408     }
0409     if (bspin == -1) {
0410     hists->phi_pT_blue_down->FillHists(pT, phiblue);
0411     if (eta > 0) hists->phi_pT_blue_down_forward->FillHists(pT, phiblue);
0412     if (eta < 0) hists->phi_pT_blue_down_backward->FillHists(pT, phiblue);
0413     hists->phi_xF_blue_down->FillHists(xFblue, phiblue);
0414     hists->phi_eta_blue_down->FillHists(eta, phiblue);
0415     hists->phi_vtxz_blue_down->FillHists(vtxz, phiblue);
0416     if (abs(eta) < 0.35)
0417     {
0418         hists_lowEta->phi_pT_blue_down->FillHists(pT, phiblue);
0419         if (eta > 0) hists_lowEta->phi_pT_blue_down_forward->FillHists(pT, phiblue);
0420         if (eta < 0) hists_lowEta->phi_pT_blue_down_backward->FillHists(pT, phiblue);
0421         hists_lowEta->phi_xF_blue_down->FillHists(xFblue, phiblue);
0422         hists_lowEta->phi_eta_blue_down->FillHists(eta, phiblue);
0423         hists_lowEta->phi_vtxz_blue_down->FillHists(vtxz, phiblue);
0424     }
0425     if (abs(eta) > 0.35)
0426     {
0427         hists_highEta->phi_pT_blue_down->FillHists(pT, phiblue);
0428         if (eta > 0) hists_highEta->phi_pT_blue_down_forward->FillHists(pT, phiblue);
0429         if (eta < 0) hists_highEta->phi_pT_blue_down_backward->FillHists(pT, phiblue);
0430         hists_highEta->phi_xF_blue_down->FillHists(xFblue, phiblue);
0431         hists_highEta->phi_eta_blue_down->FillHists(eta, phiblue);
0432         hists_highEta->phi_vtxz_blue_down->FillHists(vtxz, phiblue);
0433     }
0434     }
0435     if (yspin == 1) {
0436     hists->phi_pT_yellow_up->FillHists(pT, phiyellow);
0437     if (eta < 0) hists->phi_pT_yellow_up_forward->FillHists(pT, phiyellow);
0438     if (eta > 0) hists->phi_pT_yellow_up_backward->FillHists(pT, phiyellow);
0439     hists->phi_xF_yellow_up->FillHists(xFyellow, phiyellow);
0440     hists->phi_eta_yellow_up->FillHists(eta, phiyellow);
0441     hists->phi_vtxz_yellow_up->FillHists(vtxz, phiyellow);
0442     if (abs(eta) < 0.35)
0443     {
0444         hists_lowEta->phi_pT_yellow_up->FillHists(pT, phiyellow);
0445         if (eta < 0) hists_lowEta->phi_pT_yellow_up_forward->FillHists(pT, phiyellow);
0446         if (eta > 0) hists_lowEta->phi_pT_yellow_up_backward->FillHists(pT, phiyellow);
0447         hists_lowEta->phi_xF_yellow_up->FillHists(xFyellow, phiyellow);
0448         hists_lowEta->phi_eta_yellow_up->FillHists(eta, phiyellow);
0449         hists_lowEta->phi_vtxz_yellow_up->FillHists(vtxz, phiyellow);
0450     }
0451     if (abs(eta) > 0.35)
0452     {
0453         hists_highEta->phi_pT_yellow_up->FillHists(pT, phiyellow);
0454         if (eta < 0) hists_highEta->phi_pT_yellow_up_forward->FillHists(pT, phiyellow);
0455         if (eta > 0) hists_highEta->phi_pT_yellow_up_backward->FillHists(pT, phiyellow);
0456         hists_highEta->phi_xF_yellow_up->FillHists(xFyellow, phiyellow);
0457         hists_highEta->phi_eta_yellow_up->FillHists(eta, phiyellow);
0458         hists_highEta->phi_vtxz_yellow_up->FillHists(vtxz, phiyellow);
0459     }
0460     }
0461     if (yspin == -1) {
0462     hists->phi_pT_yellow_down->FillHists(pT, phiyellow);
0463     if (eta < 0) hists->phi_pT_yellow_down_forward->FillHists(pT, phiyellow);
0464     if (eta > 0) hists->phi_pT_yellow_down_backward->FillHists(pT, phiyellow);
0465     hists->phi_xF_yellow_down->FillHists(xFyellow, phiyellow);
0466     hists->phi_eta_yellow_down->FillHists(eta, phiyellow);
0467     hists->phi_vtxz_yellow_down->FillHists(vtxz, phiyellow);
0468     if (abs(eta) < 0.35)
0469     {
0470         hists_lowEta->phi_pT_yellow_down->FillHists(pT, phiyellow);
0471         if (eta < 0) hists_lowEta->phi_pT_yellow_down_forward->FillHists(pT, phiyellow);
0472         if (eta > 0) hists_lowEta->phi_pT_yellow_down_backward->FillHists(pT, phiyellow);
0473         hists_lowEta->phi_xF_yellow_down->FillHists(xFyellow, phiyellow);
0474         hists_lowEta->phi_eta_yellow_down->FillHists(eta, phiyellow);
0475         hists_lowEta->phi_vtxz_yellow_down->FillHists(vtxz, phiyellow);
0476     }
0477     if (abs(eta) > 0.35)
0478     {
0479         hists_highEta->phi_pT_yellow_down->FillHists(pT, phiyellow);
0480         if (eta < 0) hists_highEta->phi_pT_yellow_down_forward->FillHists(pT, phiyellow);
0481         if (eta > 0) hists_highEta->phi_pT_yellow_down_backward->FillHists(pT, phiyellow);
0482         hists_highEta->phi_xF_yellow_down->FillHists(xFyellow, phiyellow);
0483         hists_highEta->phi_eta_yellow_down->FillHists(eta, phiyellow);
0484         hists_highEta->phi_vtxz_yellow_down->FillHists(vtxz, phiyellow);
0485     }
0486     }
0487 }
0488
0489 /* void TSSAhistmaker::FillAllPhiHists(int index) */
0490 /* { */
0491 /*     FillPhiHists("pi0", i); */
0492 /*     FillPhiHists("eta", i); */
0493 /*     FillPhiHists("pi0bkgr", i); */
0494 /*     FillPhiHists("etabkgr", i); */
0495 /* } */
0496
0497 /* void TSSAhistmaker::LoopClusters() */
0498 /* { */
0499 /* } */
0500
0501 /* void TSSAhistmaker::LoopDiphotons() */
0502 /* { */
0503 /* } */
0504
0505 void TSSAhistmaker::LoopTrees()
0506 {
0507     int nevents_event = tree_event->GetEntries();
0508     int nevents_clusters = tree_clusters->GetEntries();
0509     int nevents_diphotons = tree_diphotons->GetEntries();
0510     if (nevents_clusters != nevents_diphotons) {
0511     std::cout << "Found different number of events in Cluster and Diphoton! Exiting!" << std::endl;
0512     exit(1);
0513     }
0514     if (nevents_event != nevents_diphotons) {
0515     std::cout << "Found different number of events in Event and Diphoton! Exiting!" << std::endl;
0516     exit(1);
0517     }
0518     int nevents = nevents_clusters;
0519     std::cout << "Found " << nevents << " total events" << std::endl;
0520     for (int i=0; i<nevents; i++) {
0521     if (i%10000 == 0) std::cout << "Event " << i << std::endl;
0522     h_Nevents->Fill(0);
0523     tree_event->GetEntry(i);
0524     tree_clusters->GetEntry(i);
0525
0526     int nclusters = cluster_E->size();
0527     h_vtxz->Fill(vtxz);
0528     /* h_nAllClusters->Fill(nAllClusters); */
0529     /* h_nGoodClusters->Fill(nclusters); */
0530     if (minbiastrig_scaled) h_nClustersMinbiasTrig->Fill(nclusters);
0531     if (photontrig_scaled) h_nClustersPhotonTrig->Fill(nclusters);
0532
0533     // Armenteros-Podolanski plot
0534     if (nclusters < 2) continue;
0535     for (int j=0; j<(nclusters-1); j++) {
0536         for (int k=(j+1); k<nclusters; k++) {
0537         double E1 = cluster_E->at(j);
0538         double E2 = cluster_E->at(k);
0539         double asymE = (E1 - E2)/(E1 + E2);
0540         double eta1 = cluster_Eta->at(j);
0541         double eta2 = cluster_Eta->at(k);
0542         double phi1 = cluster_Phi->at(j);
0543         double phi2 = cluster_Phi->at(k);
0544         double pT1 = cluster_pT->at(j);
0545         double pT2 = cluster_pT->at(k);
0546         TLorentzVector v1, v2;
0547         v1.SetPtEtaPhiE(pT1, eta1, phi1, E1);
0548         v2.SetPtEtaPhiE(pT2, eta2, phi2, E2);
0549         TLorentzVector diphoton = v1 + v2;
0550         double mass = diphoton.M();
0551         double pT = diphoton.Pt();
0552         if (pT < 1.0) continue;
0553         double deltaR = v1.DeltaR(v2);
0554         TVector3 v1_3d = v1.Vect();
0555         TVector3 v2_3d = v2.Vect();
0556         TVector3 d_3d = diphoton.Vect();
0557         double p_T1 = v1_3d.Pt(d_3d);
0558         double p_T2 = v2_3d.Pt(d_3d);
0559         double p_L1 = TMath::Sqrt(v1_3d.Mag2() - p_T1*p_T1);
0560         double p_L2 = TMath::Sqrt(v2_3d.Mag2() - p_T2*p_T2);
0561         double alpha = (p_L1 - p_L2)/(p_L1 + p_L2);
0562         h_armen_all->Fill(alpha, p_T1);
0563         if (InRange(mass, pi0MassRange)) h_armen_pi0->Fill(alpha, p_T1);
0564         if (InRange(mass, pi0BkgrLowRange) || InRange(mass, pi0BkgrHighRange)) h_armen_pi0bkgr->Fill(alpha, p_T1);
0565         if (InRange(mass, etaMassRange)) {
0566             h_armen_eta->Fill(alpha, p_T1);
0567             if (pT > 7.0) h_armen_eta_highpT->Fill(alpha, p_T1);
0568         }
0569         if (InRange(mass, etaBkgrLowRange) || InRange(mass, etaBkgrHighRange)) {
0570             h_armen_etabkgr->Fill(alpha, p_T1);
0571             if (pT > 7.0) h_armen_etabkgr_highpT->Fill(alpha, p_T1);
0572         }
0573         h_armen_p_L->Fill(p_L1);
0574         h_armen_p_L->Fill(p_L2);
0575         h_armen_p_T->Fill(p_T1);
0576         h_armen_alpha_alphaE->Fill(alpha, asymE);
0577         h_armen_alpha_deltaR->Fill(alpha, deltaR);
0578         if (InRange(mass, pi0MassRange)) h_armen_alpha_deltaR_pi0->Fill(alpha, deltaR);
0579         if (InRange(mass, etaMassRange)) h_armen_alpha_deltaR_eta->Fill(alpha, deltaR);
0580         h_armen_p_T_deltaR->Fill(p_T1, deltaR);
0581         h_armen_pT_p_L->Fill(pT, p_L1);
0582         h_armen_pT_p_L->Fill(pT, p_L2);
0583         h_armen_pT_p_T->Fill(pT, p_T1);
0584         }
0585     }
0586
0587
0588     tree_diphotons->GetEntry(i);
0589     int ndiphotons = diphoton_E->size();
0590     if (ndiphotons > 0) {
0591         for (int j=0; j<ndiphotons; j++) {
0592         // additional diphoton cuts here
0593         float mass = diphoton_M->at(j);
0594         float pT = diphoton_pT->at(j);
0595         float asym = diphoton_asym->at(j);
0596         float deltaR = diphoton_deltaR->at(j);
0597         if (pT < 1.0) continue;
0598         if (asym > 0.7) continue;
0599         float xF = diphoton_xF->at(j);
0600         float eta = diphoton_Eta->at(j);
0601         /* if (pT > 4.0 && deltaR > 1.0) continue; */
0602
0603         /* float pT1 = cluster_pT->at(diphoton_clus1index->at(j)); */
0604         /* float pT2 = cluster_pT->at(diphoton_clus2index->at(j)); */
0605         /* if (pT1 > 1.5 && pT2 > 1.5) { */
0606             /* h_diphotonMass->Fill(mass); */
0607             /* bhs_diphotonMass_pT->FillHists(pT, mass); */
0608         /* } */
0609         h_diphotonMass->Fill(mass);
0610         bhs_diphotonMass_pT_pi0binning->FillHists(pT, mass);
0611         bhs_diphotonMass_pT_etabinning->FillHists(pT, mass);
0612         if (InRange(mass, pi0MassRange)) bhs_pi0_pT_pT->FillHists(pT, pT);
0613         if (InRange(mass, etaMassRange)) bhs_eta_pT_pT->FillHists(pT, pT);
0614         bhs_diphotonxF_xF->FillHists(xF, xF);
0615         bhs_diphotoneta_eta->FillHists(eta, eta);
0616         bhs_diphotonvtxz_vtxz->FillHists(vtxz, vtxz);
0617
0618         h_DiphotonMassAsym->Fill(mass, asym);
0619         h_DiphotonMassdeltaR->Fill(mass, deltaR);
0620         h_DiphotondeltaRAsym->Fill(deltaR, asym);
0621         if (pT > 7.0) {
0622             h_DiphotonMassAsym_highpT->Fill(mass, asym);
0623             h_DiphotonMassdeltaR_highpT->Fill(mass, deltaR);
0624             h_DiphotondeltaRAsym_highpT->Fill(deltaR, asym);
0625             h_DiphotondeltaRAsym_highpT_smalldR->Fill(deltaR, asym);
0626         }
0627
0628         if (InRange(mass, pi0MassRange)) {
0629             h_Npi0s->Fill(0);
0630             FillPhiHists("pi0", j);
0631             h_DiphotondeltaRAsym_pi0->Fill(deltaR, asym);
0632             h_DiphotondeltaRAsym_pi0_smalldR->Fill(deltaR, asym);
0633         }
0634         if (InRange(mass, pi0BkgrLowRange) || InRange(mass, pi0BkgrHighRange)) {
0635             FillPhiHists("pi0bkgr", j);
0636         }
0637         if (InRange(mass, etaMassRange)) {
0638             h_Netas->Fill(0);
0639             FillPhiHists("eta", j);
0640             h_DiphotondeltaRAsym_eta->Fill(deltaR, asym);
0641             if (pT > 7.0) h_DiphotondeltaRAsym_eta_highpT->Fill(deltaR, asym);
0642         }
0643         if (InRange(mass, etaBkgrLowRange) || InRange(mass, etaBkgrHighRange)) {
0644             FillPhiHists("etabkgr", j);
0645             h_DiphotondeltaRAsym_etabkgr->Fill(deltaR, asym);
0646             if (pT > 7.0) h_DiphotondeltaRAsym_etabkgr_highpT->Fill(deltaR, asym);
0647         }
0648         }
0649     }
0650     }
0651 }
0652
0653 void TSSAhistmaker::Cleanup()
0654 {
0655     outfile_hists->Write();
0656     outfile_hists->Close();
0657     delete outfile_hists;
0658
0659     /* infile_trees->Close(); */
0660     /* delete infile_trees; */
0661
0662     delete pi0Hists; delete etaHists; delete pi0BkgrHists; delete etaBkgrHists;
0663     delete pi0Hists_lowEta; delete etaHists_lowEta; delete pi0BkgrHists_lowEta; delete etaBkgrHists_lowEta;
0664     delete pi0Hists_highEta; delete etaHists_highEta; delete pi0BkgrHists_highEta; delete etaBkgrHists_highEta;
0665 }
0666
0667 void TSSAhistmaker::main()
0668 {
0669     GetTrees();
0670     MakeAllHists();
0671     LoopTrees();
0672     Cleanup();
0673     std::cout << "All done!" << std::endl;
0674 }