From 0bb0c301ad86abdf4247bc3e65a72bb195f20b60 Mon Sep 17 00:00:00 2001 From: Chris Date: Wed, 6 Nov 2024 20:31:45 +0100 Subject: [PATCH 1/4] Initial commit, working DPt and DY binning --- .../DataMCComp.cpp | 185 ++++++++++++++++++ .../20241106_DzeroUPCDataMCComp/README.md | 3 + .../configs/dataMCCompare.config | 4 + .../20241106_DzeroUPCDataMCComp/makefile | 6 + 4 files changed, 198 insertions(+) create mode 100644 MainAnalysis/20241106_DzeroUPCDataMCComp/DataMCComp.cpp create mode 100644 MainAnalysis/20241106_DzeroUPCDataMCComp/README.md create mode 100644 MainAnalysis/20241106_DzeroUPCDataMCComp/configs/dataMCCompare.config create mode 100644 MainAnalysis/20241106_DzeroUPCDataMCComp/makefile diff --git a/MainAnalysis/20241106_DzeroUPCDataMCComp/DataMCComp.cpp b/MainAnalysis/20241106_DzeroUPCDataMCComp/DataMCComp.cpp new file mode 100644 index 00000000..f4c00050 --- /dev/null +++ b/MainAnalysis/20241106_DzeroUPCDataMCComp/DataMCComp.cpp @@ -0,0 +1,185 @@ +//Created C. McGinn, 2024.11.06 +//contact cfmcginn@mit.edu/cffionn@skype/cfmcginn@github + +//c and cpp dependencies +#include +#include +#include +#include + +//ROOT +#include "TEnv.h" +#include "TFile.h" +#include "TH1D.h" +#include "TTree.h" + +std::vector commaSepStrToDVect(std::string inStr) +{ + std::vector vectD = {}; + //Test the vector input + if(inStr.size() == 0){ + std::cout << __PRETTY_FUNCTION__ << " WARNING: Passed empty string, returning empty vector. check input" << std::endl; + return vectD; + } + + if(inStr.substr(inStr.size()-1, 1).find(",") == std::string::npos){ + inStr = inStr + ","; + } + + if(inStr.find(",") == std::string::npos){ + std::cout << __PRETTY_FUNCTION__ << " WARNING: Not a comma separated list, returning empty vector. check input" << std::endl; + return vectD; + } + + while(inStr.find(",") != std::string::npos){ + double value = std::stod(inStr.substr(0, inStr.find(","))); + vectD.push_back(value); + inStr.replace(0, inStr.find(",")+1, ""); + } + + return vectD; +} + +int DataMCComp(std::string inConfigName) +{ + //config handler - + TEnv* config_p = new TEnv(inConfigName.c_str()); + //Grab the input skim file + std::string skimFileName = config_p->GetValue("SKIMFILENAME", ""); + std::vector skimFileNames; + + //Process if its a root or a txt list of root files + if(skimFileName.find(".root") != std::string::npos) skimFileNames.push_back(skimFileName); + else if(skimFileName.find(".txt") != std::string::npos){ + std::ifstream inFile(skimFileName); + std::string tempStr; + while(std::getline(inFile, tempStr)){ + if(tempStr.find(".root") != std::string::npos) skimFileNames.push_back(skimFileName); + } + inFile.close(); + } + + //grab dPtBins + std::vector dPtBins = commaSepStrToDVect(config_p->GetValue("DPTBINS", "")); + if(dPtBins.size() == 0) return 1; + + //grab dYBins + std::vector dYBins = commaSepStrToDVect(config_p->GetValue("DYBINS", "")); + if(dYBins.size() == 0) return 1; + + //Grab the outputfile name, if it exist; else default name + std::string outFileName = config_p->GetValue("OUTFILENAME", "out.root"); + TFile* outFile_p = new TFile(outFileName.c_str(), "RECREATE"); + + const Int_t nMaxBins = 20; + + const Int_t nDPtBins = dPtBins.size()-1; + const Int_t nDYBins = dYBins.size()-1; + + TH1D* dPt_h[nMaxBins]; + TH1D* dY_h[nMaxBins]; + //Pt histograms in y bins + for(unsigned int yI = 0; yI < nDYBins; ++yI){ + std::string histName = Form("dPt_YBin%d_h", yI); + std::string histTitle = Form(";D^{0} p_{T} (GeV), %.1f < y < %.1f; Counts", dYBins[yI], dYBins[yI+1]); + + dPt_h[yI] = new TH1D(histName.c_str(), histTitle.c_str(), dPtBins.size()-1, dPtBins.data()); + } + std::string histName = "dPt_YBinAll_h"; + std::string histTitle = Form(";D^{0} p_{T} (GeV), %.1f < y < %.1f; Counts", dYBins[0], dYBins[nDYBins]); + dPt_h[nDYBins] = new TH1D(histName.c_str(), histTitle.c_str(), dPtBins.size()-1, dPtBins.data()); + + //Y histograms in pt bins + for(unsigned int pI = 0; pI < nDPtBins; ++pI){ + histName = Form("dY_PtBin%d_h", pI); + histTitle = Form(";D^{0} y, %.1f < p_{T} < %.1f (GeV); Counts", dPtBins[pI], dPtBins[pI+1]); + + dY_h[pI] = new TH1D(histName.c_str(), histTitle.c_str(), dYBins.size()-1, dYBins.data()); + } + histName = "dY_PtBinAll_h"; + histTitle = Form(";D^{0} y, %.1f < p_{T} < %.1f (GeV); Counts", dPtBins[0], dPtBins[nDPtBins]); + dY_h[nDPtBins] = new TH1D(histName.c_str(), histTitle.c_str(), dYBins.size()-1, dYBins.data()); + + //Done handling config input + + //Process inputs + std::cout << "Processing " << skimFileNames.size() << " skim files...." << std::endl; + for(unsigned int fI = 0; fI < skimFileNames.size(); ++fI){ + std::cout << " File " << fI << "/" << skimFileNames.size() << " (" << ((Double_t)fI)/(Double_t)skimFileNames.size() << "): " << skimFileNames[fI] << std::endl; + + //Define the variables to be used w/ skim tree + std::vector* Dpt_p=nullptr; + std::vector* Dy_p=nullptr; + + //Grab the input file and Tree + TFile* inFile_p = new TFile(skimFileNames[fI].c_str(), "READ"); + TTree* Tree = (TTree*)inFile_p->Get("Tree"); + + //Tree branch handling + Tree->SetBranchStatus("*", 0); + Tree->SetBranchStatus("Dpt", 1); + Tree->SetBranchStatus("Dy", 1); + + Tree->SetBranchAddress("Dpt", &Dpt_p); + Tree->SetBranchAddress("Dy", &Dy_p); + + //Loop over entries + ULong64_t nEntries = Tree->GetEntries(); + for(ULong64_t entry = 0; entry < nEntries; ++entry){ + Tree->GetEntry(entry); + for(unsigned int dI = 0; dI < Dpt_p->size(); ++dI){ + + int yPos = -1; + for(Int_t yI = 0; yI < nDYBins; ++yI){ + if(Dy_p->at(dI) >= dYBins[yI] && Dy_p->at(dI) <= dYBins[yI+1]){ + yPos = yI; + break; + } + } + + if(yPos >= 0){ + dPt_h[yPos]->Fill(Dpt_p->at(dI)); + dPt_h[nDYBins]->Fill(Dpt_p->at(dI)); + } + } + } + + inFile_p->Close(); + delete inFile_p; + }//complete file loop + + //outfile handling + outFile_p->cd(); + + for(unsigned int yI = 0; yI < nDYBins+1; ++yI){ + dPt_h[yI]->Write("", TObject::kOverwrite); + delete dPt_h[yI]; + } + + for(unsigned int pI = 0; pI < nDPtBins; ++pI){ + dY_h[pI]->Write("", TObject::kOverwrite); + delete dY_h[pI]; + } + + //Write out the config so there is a record of production + config_p->Write("config", TObject::kOverwrite); + delete config_p; + + outFile_p->Close(); + delete outFile_p; + + std::cout << "DataMCComp complete! return 0" << std::endl; + return 0; +} + +int main(int argc, char* argv[]) +{ + if(argc != 2){ + std::cout << "Usage: ./DataMCComp . return 1" << std::endl; + return 1; + } + + int retVal = 0; + retVal += DataMCComp(argv[1]); + return retVal; +} diff --git a/MainAnalysis/20241106_DzeroUPCDataMCComp/README.md b/MainAnalysis/20241106_DzeroUPCDataMCComp/README.md new file mode 100644 index 00000000..8bac8390 --- /dev/null +++ b/MainAnalysis/20241106_DzeroUPCDataMCComp/README.md @@ -0,0 +1,3 @@ +To build, do + +make \ No newline at end of file diff --git a/MainAnalysis/20241106_DzeroUPCDataMCComp/configs/dataMCCompare.config b/MainAnalysis/20241106_DzeroUPCDataMCComp/configs/dataMCCompare.config new file mode 100644 index 00000000..ab4870be --- /dev/null +++ b/MainAnalysis/20241106_DzeroUPCDataMCComp/configs/dataMCCompare.config @@ -0,0 +1,4 @@ +SKIMFILENAME: /home/data/public/yuchenc/SkimsCheck_20241104/NewSkim/exampleData.root + +DPTBINS: 0,2,5,8,12,20 +DYBINS: -2,-1,0,1,2 \ No newline at end of file diff --git a/MainAnalysis/20241106_DzeroUPCDataMCComp/makefile b/MainAnalysis/20241106_DzeroUPCDataMCComp/makefile new file mode 100644 index 00000000..11b70565 --- /dev/null +++ b/MainAnalysis/20241106_DzeroUPCDataMCComp/makefile @@ -0,0 +1,6 @@ +default: DataMCComp + +DataMCComp: DataMCComp.cpp + g++ DataMCComp.cpp -O3 -I. -I$(ProjectBase)/CommonCode/include -I./include \ + -o DataMCComp `root-config --cflags --glibs` \ + $(ProjectBase)/CommonCode/library/Messenger.o From 32d341ae3602238e7e21d24d5f93b733c4caafca Mon Sep 17 00:00:00 2001 From: Chris Date: Thu, 7 Nov 2024 11:37:31 +0100 Subject: [PATCH 2/4] Updated for some bin protections, clarifying this is raw D, etc --- .../DataMCComp.cpp | 116 +++++++++++------- 1 file changed, 74 insertions(+), 42 deletions(-) diff --git a/MainAnalysis/20241106_DzeroUPCDataMCComp/DataMCComp.cpp b/MainAnalysis/20241106_DzeroUPCDataMCComp/DataMCComp.cpp index f4c00050..ff00c102 100644 --- a/MainAnalysis/20241106_DzeroUPCDataMCComp/DataMCComp.cpp +++ b/MainAnalysis/20241106_DzeroUPCDataMCComp/DataMCComp.cpp @@ -13,6 +13,9 @@ #include "TH1D.h" #include "TTree.h" +//Local CommonCode +#include "Messenger.h" + std::vector commaSepStrToDVect(std::string inStr) { std::vector vectD = {}; @@ -40,6 +43,34 @@ std::vector commaSepStrToDVect(std::string inStr) return vectD; } +int getPosFromBins(double val, std::vector bins) +{ + int retPos = -1; + if(bins.size() <= 1){ + std::cout << __PRETTY_FUNCTION__ << " WARNING: Given bins.size() = " << bins.size() << ", return -1" << std::endl; + return retPos; + } + + for(unsigned int bI = 0; bI < bins.size()-1; ++bI){ + if(val >= bins[bI] && val < bins[bI+1]){ + retPos = bI; + break; + } + } + + return retPos; +} + +bool binsMaxCheck(std::string binStr, int nBins, int nMaxBins) +{ + if(nBins > nMaxBins){ + std::cout << __PRETTY_FUNCTION__ << " WARNING: Given nBins=" << nBins << " for '' exceeds max. return false" << std::endl; + return false; + } + + return true; +} + int DataMCComp(std::string inConfigName) { //config handler - @@ -72,33 +103,39 @@ int DataMCComp(std::string inConfigName) TFile* outFile_p = new TFile(outFileName.c_str(), "RECREATE"); const Int_t nMaxBins = 20; - const Int_t nDPtBins = dPtBins.size()-1; const Int_t nDYBins = dYBins.size()-1; + + if(!binsMaxCheck("DPt", nDPtBins, nMaxBins)) return 1; + if(!binsMaxCheck("DY", nDYBins, nMaxBins)) return 1; + - TH1D* dPt_h[nMaxBins]; - TH1D* dY_h[nMaxBins]; + TH1D* rawDPt_h[nMaxBins]; + TH1D* rawDY_h[nMaxBins]; + TH1D* rawDM_h[nMaxBins][nMaxBins]; + //Pt histograms in y bins for(unsigned int yI = 0; yI < nDYBins; ++yI){ - std::string histName = Form("dPt_YBin%d_h", yI); - std::string histTitle = Form(";D^{0} p_{T} (GeV), %.1f < y < %.1f; Counts", dYBins[yI], dYBins[yI+1]); + std::string histName = Form("rawDPt_YBin%d_h", yI); + std::string histTitle = Form(";D^{0} cand. p_{T} (GeV), %.1f < y < %.1f; Counts", dYBins[yI], dYBins[yI+1]); - dPt_h[yI] = new TH1D(histName.c_str(), histTitle.c_str(), dPtBins.size()-1, dPtBins.data()); + rawDPt_h[yI] = new TH1D(histName.c_str(), histTitle.c_str(), dPtBins.size()-1, dPtBins.data()); } - std::string histName = "dPt_YBinAll_h"; - std::string histTitle = Form(";D^{0} p_{T} (GeV), %.1f < y < %.1f; Counts", dYBins[0], dYBins[nDYBins]); - dPt_h[nDYBins] = new TH1D(histName.c_str(), histTitle.c_str(), dPtBins.size()-1, dPtBins.data()); + std::string histName = "rawDPt_YBinAll_h"; + std::string histTitle = Form(";D^{0} cand. p_{T} (GeV), %.1f < y < %.1f; Counts", dYBins[0], dYBins[nDYBins]); + rawDPt_h[nDYBins] = new TH1D(histName.c_str(), histTitle.c_str(), dPtBins.size()-1, dPtBins.data()); //Y histograms in pt bins for(unsigned int pI = 0; pI < nDPtBins; ++pI){ - histName = Form("dY_PtBin%d_h", pI); - histTitle = Form(";D^{0} y, %.1f < p_{T} < %.1f (GeV); Counts", dPtBins[pI], dPtBins[pI+1]); + histName = Form("rawDY_PtBin%d_h", pI); + histTitle = Form(";D^{0} cand. y, %.1f < p_{T} < %.1f (GeV); Counts", dPtBins[pI], dPtBins[pI+1]); - dY_h[pI] = new TH1D(histName.c_str(), histTitle.c_str(), dYBins.size()-1, dYBins.data()); + rawDY_h[pI] = new TH1D(histName.c_str(), histTitle.c_str(), dYBins.size()-1, dYBins.data()); } - histName = "dY_PtBinAll_h"; - histTitle = Form(";D^{0} y, %.1f < p_{T} < %.1f (GeV); Counts", dPtBins[0], dPtBins[nDPtBins]); - dY_h[nDPtBins] = new TH1D(histName.c_str(), histTitle.c_str(), dYBins.size()-1, dYBins.data()); + histName = "rawDY_PtBinAll_h"; + histTitle = Form(";D^{0} cand. y, %.1f < p_{T} < %.1f (GeV); Counts", dPtBins[0], dPtBins[nDPtBins]); + rawDY_h[nDPtBins] = new TH1D(histName.c_str(), histTitle.c_str(), dYBins.size()-1, dYBins.data()); + //Done handling config input @@ -113,36 +150,31 @@ int DataMCComp(std::string inConfigName) //Grab the input file and Tree TFile* inFile_p = new TFile(skimFileNames[fI].c_str(), "READ"); - TTree* Tree = (TTree*)inFile_p->Get("Tree"); - - //Tree branch handling - Tree->SetBranchStatus("*", 0); - Tree->SetBranchStatus("Dpt", 1); - Tree->SetBranchStatus("Dy", 1); - - Tree->SetBranchAddress("Dpt", &Dpt_p); - Tree->SetBranchAddress("Dy", &Dy_p); + //Use UPC DZeroTree + DzeroUPCTreeMessenger *MDzeroUPC = new DzeroUPCTreeMessenger(inFile_p, "Tree"); //Loop over entries - ULong64_t nEntries = Tree->GetEntries(); + ULong64_t nEntries = MDzeroUPC->GetEntries(); for(ULong64_t entry = 0; entry < nEntries; ++entry){ - Tree->GetEntry(entry); - for(unsigned int dI = 0; dI < Dpt_p->size(); ++dI){ - - int yPos = -1; - for(Int_t yI = 0; yI < nDYBins; ++yI){ - if(Dy_p->at(dI) >= dYBins[yI] && Dy_p->at(dI) <= dYBins[yI+1]){ - yPos = yI; - break; - } - } + MDzeroUPC->GetEntry(entry); + for(unsigned int dI = 0; dI < MDzeroUPC->Dpt->size(); ++dI){ + int yPos = getPosFromBins(MDzeroUPC->Dy->at(dI), dYBins); + int ptPos = getPosFromBins(MDzeroUPC->Dpt->at(dI), dPtBins); + if(yPos >= 0){ - dPt_h[yPos]->Fill(Dpt_p->at(dI)); - dPt_h[nDYBins]->Fill(Dpt_p->at(dI)); + rawDPt_h[yPos]->Fill(MDzeroUPC->Dpt->at(dI)); + rawDPt_h[nDYBins]->Fill(MDzeroUPC->Dpt->at(dI)); + } + + if(ptPos >= 0){ + rawDY_h[ptPos]->Fill(MDzeroUPC->Dy->at(dI)); + rawDY_h[nDPtBins]->Fill(MDzeroUPC->Dy->at(dI)); } } } + + delete MDzeroUPC; inFile_p->Close(); delete inFile_p; @@ -152,13 +184,13 @@ int DataMCComp(std::string inConfigName) outFile_p->cd(); for(unsigned int yI = 0; yI < nDYBins+1; ++yI){ - dPt_h[yI]->Write("", TObject::kOverwrite); - delete dPt_h[yI]; + rawDPt_h[yI]->Write("", TObject::kOverwrite); + delete rawDPt_h[yI]; } - for(unsigned int pI = 0; pI < nDPtBins; ++pI){ - dY_h[pI]->Write("", TObject::kOverwrite); - delete dY_h[pI]; + for(unsigned int pI = 0; pI < nDPtBins+1; ++pI){ + rawDY_h[pI]->Write("", TObject::kOverwrite); + delete rawDY_h[pI]; } //Write out the config so there is a record of production From a9b2e19f8b2aff36942473be1964d853a8d149a7 Mon Sep 17 00:00:00 2001 From: Chris Date: Thu, 7 Nov 2024 20:49:13 +0100 Subject: [PATCH 3/4] Added processing timers and DMass --- .../DataMCComp.cpp | 115 +++++++++++++----- .../configs/dataMCCompare.config | 10 +- 2 files changed, 94 insertions(+), 31 deletions(-) diff --git a/MainAnalysis/20241106_DzeroUPCDataMCComp/DataMCComp.cpp b/MainAnalysis/20241106_DzeroUPCDataMCComp/DataMCComp.cpp index ff00c102..d649262d 100644 --- a/MainAnalysis/20241106_DzeroUPCDataMCComp/DataMCComp.cpp +++ b/MainAnalysis/20241106_DzeroUPCDataMCComp/DataMCComp.cpp @@ -71,6 +71,34 @@ bool binsMaxCheck(std::string binStr, int nBins, int nMaxBins) return true; } +void bookHist1D(std::string nameStr, std::string titleStr, std::vector histBins, std::vector axisBins, TH1D* hist[]) +{ + for(unsigned int bI = 0; bI < histBins.size()-1; ++bI){ + std::string localName = nameStr + std::to_string(bI) + "_h"; + std::string localTitle = Form(titleStr.c_str(), histBins[bI], histBins[bI+1]); + + hist[bI] = new TH1D(localName.c_str(), localTitle.c_str(), axisBins.size()-1, axisBins.data()); + } + std::string localName = nameStr + "All_h"; + std::string localTitle = Form(titleStr.c_str(), histBins[0], histBins[histBins.size()-1]); + hist[histBins.size()-1] = new TH1D(localName.c_str(), localTitle.c_str(), axisBins.size()-1, axisBins.data()); + + return; +} + +std::vector getLinBins(int nbins, double low, double high) +{ + std::vector vect; + double binWidth = (high - low)/(double)nbins; + vect.push_back(low); + for(int bI = 1; bI < nbins; ++bI){ + vect.push_back(low+binWidth*(double)bI); + } + vect.push_back(high); + + return vect; +} + int DataMCComp(std::string inConfigName) { //config handler - @@ -89,7 +117,7 @@ int DataMCComp(std::string inConfigName) } inFile.close(); } - + //grab dPtBins std::vector dPtBins = commaSepStrToDVect(config_p->GetValue("DPTBINS", "")); if(dPtBins.size() == 0) return 1; @@ -97,6 +125,12 @@ int DataMCComp(std::string inConfigName) //grab dYBins std::vector dYBins = commaSepStrToDVect(config_p->GetValue("DYBINS", "")); if(dYBins.size() == 0) return 1; + + //Grab dMBins + Int_t nDMBins = config_p->GetValue("NDMBINS", -1);//nonsense defaults for auto-failure + Double_t dMBinsLow = config_p->GetValue("DMBINSLOW", -1.0); + Double_t dMBinsHigh = config_p->GetValue("DMBINSHIGH", -9999.0); + std::vector dMBins = getLinBins(nDMBins, dMBinsLow, dMBinsHigh); //Grab the outputfile name, if it exist; else default name std::string outFileName = config_p->GetValue("OUTFILENAME", "out.root"); @@ -108,46 +142,52 @@ int DataMCComp(std::string inConfigName) if(!binsMaxCheck("DPt", nDPtBins, nMaxBins)) return 1; if(!binsMaxCheck("DY", nDYBins, nMaxBins)) return 1; - TH1D* rawDPt_h[nMaxBins]; TH1D* rawDY_h[nMaxBins]; TH1D* rawDM_h[nMaxBins][nMaxBins]; - //Pt histograms in y bins - for(unsigned int yI = 0; yI < nDYBins; ++yI){ - std::string histName = Form("rawDPt_YBin%d_h", yI); - std::string histTitle = Form(";D^{0} cand. p_{T} (GeV), %.1f < y < %.1f; Counts", dYBins[yI], dYBins[yI+1]); - - rawDPt_h[yI] = new TH1D(histName.c_str(), histTitle.c_str(), dPtBins.size()-1, dPtBins.data()); - } - std::string histName = "rawDPt_YBinAll_h"; - std::string histTitle = Form(";D^{0} cand. p_{T} (GeV), %.1f < y < %.1f; Counts", dYBins[0], dYBins[nDYBins]); - rawDPt_h[nDYBins] = new TH1D(histName.c_str(), histTitle.c_str(), dPtBins.size()-1, dPtBins.data()); - - //Y histograms in pt bins - for(unsigned int pI = 0; pI < nDPtBins; ++pI){ - histName = Form("rawDY_PtBin%d_h", pI); - histTitle = Form(";D^{0} cand. y, %.1f < p_{T} < %.1f (GeV); Counts", dPtBins[pI], dPtBins[pI+1]); - - rawDY_h[pI] = new TH1D(histName.c_str(), histTitle.c_str(), dYBins.size()-1, dYBins.data()); + bookHist1D("rawDPt_YBin", ";D^{0} cand. p_{T} (GeV), %.1f < y < %.1f; Counts", dYBins, dPtBins, rawDPt_h); + bookHist1D("rawDY_PtBin", ";D^{0} cand. y, %.1f < p_{T} < %.1f (GeV); Counts", dPtBins, dYBins, rawDY_h); + + //For mass, double binned, so iterate over Pt, then Y + for(Int_t pI = 0; pI < nDPtBins+1; ++pI){ + std::string histName = "rawDM_PtBinAll"; + std::string histTitle = Form(";D^{0} cand. Mass (GeV), %.1f < p_{T} < %.1f (GeV),", dPtBins[0], dPtBins[dPtBins.size()-1]); + if(pI < nDPtBins){ + histName = "rawDM_PtBin" + std::to_string(pI); + histTitle = Form(";D^{0} cand. Mass (GeV), %.1f < p_{T} < %.1f (GeV),", dPtBins[pI], dPtBins[pI+1]); + } + + histName = histName + "_Y"; + histTitle = histTitle + " %.1f < y < %.1f;Counts"; + bookHist1D(histName, histTitle, dYBins, dMBins, rawDM_h[pI]); } - histName = "rawDY_PtBinAll_h"; - histTitle = Form(";D^{0} cand. y, %.1f < p_{T} < %.1f (GeV); Counts", dPtBins[0], dPtBins[nDPtBins]); - rawDY_h[nDPtBins] = new TH1D(histName.c_str(), histTitle.c_str(), dYBins.size()-1, dYBins.data()); + const ULong64_t nEntriesToPrint = TMath::Max(1, config_p->GetValue("NENTRIESTOPRINT", 10000)); //Done handling config input + ULong64_t totalNEntries = 0; + + std::cout << "Pre-processing " << skimFileNames.size() << " skim files...." << std::endl; + for(unsigned int fI = 0; fI < skimFileNames.size(); ++fI){ + TFile* inFile_p = new TFile(skimFileNames[fI].c_str(), "READ"); + DzeroUPCTreeMessenger *MDzeroUPC = new DzeroUPCTreeMessenger(inFile_p, "Tree"); + totalNEntries += MDzeroUPC->GetEntries(); + delete MDzeroUPC; + inFile_p->Close(); + delete inFile_p; + } + //Process inputs std::cout << "Processing " << skimFileNames.size() << " skim files...." << std::endl; - for(unsigned int fI = 0; fI < skimFileNames.size(); ++fI){ - std::cout << " File " << fI << "/" << skimFileNames.size() << " (" << ((Double_t)fI)/(Double_t)skimFileNames.size() << "): " << skimFileNames[fI] << std::endl; + std::cout << "Total entries: " << totalNEntries << "..." << std::endl; - //Define the variables to be used w/ skim tree - std::vector* Dpt_p=nullptr; - std::vector* Dy_p=nullptr; - + ULong64_t nCurrEntries = 0; + + + for(unsigned int fI = 0; fI < skimFileNames.size(); ++fI){ //Grab the input file and Tree TFile* inFile_p = new TFile(skimFileNames[fI].c_str(), "READ"); //Use UPC DZeroTree @@ -156,6 +196,9 @@ int DataMCComp(std::string inConfigName) //Loop over entries ULong64_t nEntries = MDzeroUPC->GetEntries(); for(ULong64_t entry = 0; entry < nEntries; ++entry){ + if(nCurrEntries%nEntriesToPrint == 0) std::cout << " Entry " << nCurrEntries << "/" << totalNEntries << " = " << Form("%.1f", ((Double_t)100*nCurrEntries)/((Double_t)totalNEntries)) << "%, File " << fI << "/" << skimFileNames.size() << std::endl; + ++nCurrEntries; + MDzeroUPC->GetEntry(entry); for(unsigned int dI = 0; dI < MDzeroUPC->Dpt->size(); ++dI){ @@ -171,6 +214,13 @@ int DataMCComp(std::string inConfigName) rawDY_h[ptPos]->Fill(MDzeroUPC->Dy->at(dI)); rawDY_h[nDPtBins]->Fill(MDzeroUPC->Dy->at(dI)); } + + if(ptPos >= 0 && yPos >= 0){ + rawDM_h[ptPos][yPos]->Fill(MDzeroUPC->Dmass->at(dI)); + rawDM_h[nDPtBins][yPos]->Fill(MDzeroUPC->Dmass->at(dI)); + rawDM_h[ptPos][nDYBins]->Fill(MDzeroUPC->Dmass->at(dI)); + rawDM_h[nDPtBins][nDYBins]->Fill(MDzeroUPC->Dmass->at(dI)); + } } } @@ -192,7 +242,14 @@ int DataMCComp(std::string inConfigName) rawDY_h[pI]->Write("", TObject::kOverwrite); delete rawDY_h[pI]; } - + + for(Int_t pI = 0; pI < nDPtBins+1; ++pI){ + for(unsigned int yI = 0; yI < nDYBins+1; ++yI){ + rawDM_h[pI][yI]->Write("", TObject::kOverwrite); + delete rawDM_h[pI][yI]; + } + } + //Write out the config so there is a record of production config_p->Write("config", TObject::kOverwrite); delete config_p; diff --git a/MainAnalysis/20241106_DzeroUPCDataMCComp/configs/dataMCCompare.config b/MainAnalysis/20241106_DzeroUPCDataMCComp/configs/dataMCCompare.config index ab4870be..120e7065 100644 --- a/MainAnalysis/20241106_DzeroUPCDataMCComp/configs/dataMCCompare.config +++ b/MainAnalysis/20241106_DzeroUPCDataMCComp/configs/dataMCCompare.config @@ -1,4 +1,10 @@ -SKIMFILENAME: /home/data/public/yuchenc/SkimsCheck_20241104/NewSkim/exampleData.root +SKIMFILENAME: /data/NewSkims23_24/20241106_SkimOldReco23sample_DataAll.root DPTBINS: 0,2,5,8,12,20 -DYBINS: -2,-1,0,1,2 \ No newline at end of file +DYBINS: -2,-1,0,1,2 + +NDMBINS: 50 +DMBINSLOW: 1.6 +DMBINSHIGH: 2.1 + +NENTRIESTOPRINT: 1000000 \ No newline at end of file From 3682deee7bb87ff32296bb98d8ca8800b5d4c91c Mon Sep 17 00:00:00 2001 From: Chris Date: Fri, 8 Nov 2024 11:41:39 +0100 Subject: [PATCH 4/4] Updated for GammaN/NGamma binning, etc. --- .../DataMCComp.cpp | 196 ++++++++++++------ .../configs/dataMCCompare.config | 5 +- 2 files changed, 139 insertions(+), 62 deletions(-) diff --git a/MainAnalysis/20241106_DzeroUPCDataMCComp/DataMCComp.cpp b/MainAnalysis/20241106_DzeroUPCDataMCComp/DataMCComp.cpp index d649262d..6e72e7a2 100644 --- a/MainAnalysis/20241106_DzeroUPCDataMCComp/DataMCComp.cpp +++ b/MainAnalysis/20241106_DzeroUPCDataMCComp/DataMCComp.cpp @@ -28,7 +28,7 @@ std::vector commaSepStrToDVect(std::string inStr) if(inStr.substr(inStr.size()-1, 1).find(",") == std::string::npos){ inStr = inStr + ","; } - + if(inStr.find(",") == std::string::npos){ std::cout << __PRETTY_FUNCTION__ << " WARNING: Not a comma separated list, returning empty vector. check input" << std::endl; return vectD; @@ -39,7 +39,7 @@ std::vector commaSepStrToDVect(std::string inStr) vectD.push_back(value); inStr.replace(0, inStr.find(",")+1, ""); } - + return vectD; } @@ -50,7 +50,7 @@ int getPosFromBins(double val, std::vector bins) std::cout << __PRETTY_FUNCTION__ << " WARNING: Given bins.size() = " << bins.size() << ", return -1" << std::endl; return retPos; } - + for(unsigned int bI = 0; bI < bins.size()-1; ++bI){ if(val >= bins[bI] && val < bins[bI+1]){ retPos = bI; @@ -95,13 +95,13 @@ std::vector getLinBins(int nbins, double low, double high) vect.push_back(low+binWidth*(double)bI); } vect.push_back(high); - + return vect; } int DataMCComp(std::string inConfigName) { - //config handler - + //config handler - TEnv* config_p = new TEnv(inConfigName.c_str()); //Grab the input skim file std::string skimFileName = config_p->GetValue("SKIMFILENAME", ""); @@ -117,7 +117,7 @@ int DataMCComp(std::string inConfigName) } inFile.close(); } - + //grab dPtBins std::vector dPtBins = commaSepStrToDVect(config_p->GetValue("DPTBINS", "")); if(dPtBins.size() == 0) return 1; @@ -131,45 +131,84 @@ int DataMCComp(std::string inConfigName) Double_t dMBinsLow = config_p->GetValue("DMBINSLOW", -1.0); Double_t dMBinsHigh = config_p->GetValue("DMBINSHIGH", -9999.0); std::vector dMBins = getLinBins(nDMBins, dMBinsLow, dMBinsHigh); - + //Grab the outputfile name, if it exist; else default name std::string outFileName = config_p->GetValue("OUTFILENAME", "out.root"); TFile* outFile_p = new TFile(outFileName.c_str(), "RECREATE"); + //Bin handling const Int_t nMaxBins = 20; + const Int_t nGammaBins = 2; const Int_t nDPtBins = dPtBins.size()-1; const Int_t nDYBins = dYBins.size()-1; + //bin check, return if failed if(!binsMaxCheck("DPt", nDPtBins, nMaxBins)) return 1; if(!binsMaxCheck("DY", nDYBins, nMaxBins)) return 1; - - TH1D* rawDPt_h[nMaxBins]; - TH1D* rawDY_h[nMaxBins]; - TH1D* rawDM_h[nMaxBins][nMaxBins]; - - bookHist1D("rawDPt_YBin", ";D^{0} cand. p_{T} (GeV), %.1f < y < %.1f; Counts", dYBins, dPtBins, rawDPt_h); - bookHist1D("rawDY_PtBin", ";D^{0} cand. y, %.1f < p_{T} < %.1f (GeV); Counts", dPtBins, dYBins, rawDY_h); - - //For mass, double binned, so iterate over Pt, then Y - for(Int_t pI = 0; pI < nDPtBins+1; ++pI){ - std::string histName = "rawDM_PtBinAll"; - std::string histTitle = Form(";D^{0} cand. Mass (GeV), %.1f < p_{T} < %.1f (GeV),", dPtBins[0], dPtBins[dPtBins.size()-1]); - if(pI < nDPtBins){ - histName = "rawDM_PtBin" + std::to_string(pI); - histTitle = Form(";D^{0} cand. Mass (GeV), %.1f < p_{T} < %.1f (GeV),", dPtBins[pI], dPtBins[pI+1]); - } - histName = histName + "_Y"; - histTitle = histTitle + " %.1f < y < %.1f;Counts"; - bookHist1D(histName, histTitle, dYBins, dMBins, rawDM_h[pI]); + //Define gamma n and n gamma bins + std::map nGammaBinsStrMap; + nGammaBinsStrMap["NGamma"] = 0; + nGammaBinsStrMap["GammaN"] = 1; + + std::map nGammaBinsLabelMap; + nGammaBinsLabelMap["NGamma"] = "N+#gamma"; + nGammaBinsLabelMap["GammaN"] = "#gamma+N"; + + //Add info to config file + config_p->SetValue("NGAMMABINS", nGammaBins); + for(auto const& nGammaBinStr : nGammaBinsStrMap){ + std::string keyStr = "GAMMABIN." + std::to_string(nGammaBinStr.second); + std::string valStr = nGammaBinStr.first; + config_p->SetValue(keyStr.c_str(), valStr.c_str()); + + keyStr = "GAMMALABEL." + std::to_string(nGammaBinStr.second); + valStr = nGammaBinsLabelMap[nGammaBinStr.first]; + config_p->SetValue(keyStr.c_str(), valStr.c_str()); + } + + //CutPass handler + TH1D* cutPass_h = new TH1D("cutPass_h", ";Cuts;Events", 3, -0.5, 2.5); + cutPass_h->GetXaxis()->SetBinLabel(1, "No Cut"); + cutPass_h->GetXaxis()->SetBinLabel(2, "ZDC XOR"); + cutPass_h->GetXaxis()->SetBinLabel(3, "Gap Requirement"); + + //D kinematics + TH1D* rawDPt_h[nGammaBins][nMaxBins]; + TH1D* rawDY_h[nGammaBins][nMaxBins]; + TH1D* rawDM_h[nGammaBins][nMaxBins][nMaxBins]; + + for(auto const& nGammaStr: nGammaBinsStrMap){ + std::string histName = "rawDPt_" + nGammaStr.first + "_YBin"; + std::string histTitle = ";D^{0} cand. p_{T} (GeV), %.1f < y < %.1f; Counts (" + nGammaBinsLabelMap[nGammaStr.first] + ")"; + bookHist1D(histName.c_str(), histTitle.c_str(), dYBins, dPtBins, rawDPt_h[nGammaStr.second]); + + histName = "rawDY_" + nGammaStr.first + "_PtBin"; + histTitle = ";D^{0} cand. y, %.1f < p_{T} < %.1f (GeV); Counts (" + nGammaBinsLabelMap[nGammaStr.first] + ")"; + bookHist1D(histName.c_str(), histTitle.c_str(), dPtBins, dYBins, rawDY_h[nGammaStr.second]); + + //For mass, double binned, so iterate over Pt, then Y + for(Int_t pI = 0; pI < nDPtBins+1; ++pI){ + histName = "rawDM_" + nGammaStr.first + "_PtBinAll"; + histTitle = Form(";D^{0} cand. Mass (GeV), %.1f < p_{T} < %.1f (GeV),", dPtBins[0], dPtBins[dPtBins.size()-1]); + if(pI < nDPtBins){ + histName = "rawDM_" + nGammaStr.first + "_PtBin" + std::to_string(pI); + histTitle = Form(";D^{0} cand. Mass (GeV), %.1f < p_{T} < %.1f (GeV),", dPtBins[pI], dPtBins[pI+1]); + } + + histName = histName + "_Y"; + histTitle = histTitle + " %.1f < y < %.1f;Counts (" + nGammaBinsLabelMap[nGammaStr.first] + ")"; + bookHist1D(histName, histTitle, dYBins, dMBins, rawDM_h[nGammaStr.second][pI]); + } } const ULong64_t nEntriesToPrint = TMath::Max(1, config_p->GetValue("NENTRIESTOPRINT", 10000)); - + const ULong64_t nEntriesToProcess = config_p->GetValue("NENTRIESTOPROCESS", -1); + //Done handling config input ULong64_t totalNEntries = 0; - + std::cout << "Pre-processing " << skimFileNames.size() << " skim files...." << std::endl; for(unsigned int fI = 0; fI < skimFileNames.size(); ++fI){ TFile* inFile_p = new TFile(skimFileNames[fI].c_str(), "READ"); @@ -178,15 +217,21 @@ int DataMCComp(std::string inConfigName) delete MDzeroUPC; inFile_p->Close(); delete inFile_p; - } - + } + + //Check if the nentries override is in place + if(nEntriesToProcess >= 0){ + if(nEntriesToProcess >= totalNEntries){ + std::cout << "Requested nEntriesToProcess=" << nEntriesToProcess << " exceeds total number, " << totalNEntries << ". Will process total instead." << std::endl; + } + else totalNEntries = nEntriesToProcess; + } + //Process inputs std::cout << "Processing " << skimFileNames.size() << " skim files...." << std::endl; - std::cout << "Total entries: " << totalNEntries << "..." << std::endl; + std::cout << "Total entries to be processed: " << totalNEntries << "..." << std::endl; ULong64_t nCurrEntries = 0; - - for(unsigned int fI = 0; fI < skimFileNames.size(); ++fI){ //Grab the input file and Tree TFile* inFile_p = new TFile(skimFileNames[fI].c_str(), "READ"); @@ -195,68 +240,99 @@ int DataMCComp(std::string inConfigName) //Loop over entries ULong64_t nEntries = MDzeroUPC->GetEntries(); - for(ULong64_t entry = 0; entry < nEntries; ++entry){ + for(ULong64_t entry = 0; entry < nEntries; ++entry){ if(nCurrEntries%nEntriesToPrint == 0) std::cout << " Entry " << nCurrEntries << "/" << totalNEntries << " = " << Form("%.1f", ((Double_t)100*nCurrEntries)/((Double_t)totalNEntries)) << "%, File " << fI << "/" << skimFileNames.size() << std::endl; ++nCurrEntries; - + + if(nCurrEntries > totalNEntries) break; + MDzeroUPC->GetEntry(entry); + cutPass_h->Fill(0); + + //event selection handling + bool zdcGammaN = MDzeroUPC->ZDCgammaN; + bool zdcNGamma = MDzeroUPC->ZDCNgamma; + + //Impose XOR + if(zdcGammaN && zdcNGamma) continue; + if(!zdcGammaN && !zdcNGamma) continue; + + cutPass_h->Fill(1); + + //Impose gap requirement + if(zdcGammaN && !MDzeroUPC->gapgammaN) continue; + if(zdcNGamma && !MDzeroUPC->gapNgamma) continue; + + cutPass_h->Fill(2); + + std::string gammaNStr = zdcGammaN ? "GammaN" : "NGamma"; + int gammaNPos = nGammaBinsStrMap[gammaNStr]; + + //D handling for(unsigned int dI = 0; dI < MDzeroUPC->Dpt->size(); ++dI){ int yPos = getPosFromBins(MDzeroUPC->Dy->at(dI), dYBins); int ptPos = getPosFromBins(MDzeroUPC->Dpt->at(dI), dPtBins); - + if(yPos >= 0){ - rawDPt_h[yPos]->Fill(MDzeroUPC->Dpt->at(dI)); - rawDPt_h[nDYBins]->Fill(MDzeroUPC->Dpt->at(dI)); + rawDPt_h[gammaNPos][yPos]->Fill(MDzeroUPC->Dpt->at(dI)); + rawDPt_h[gammaNPos][nDYBins]->Fill(MDzeroUPC->Dpt->at(dI)); } if(ptPos >= 0){ - rawDY_h[ptPos]->Fill(MDzeroUPC->Dy->at(dI)); - rawDY_h[nDPtBins]->Fill(MDzeroUPC->Dy->at(dI)); + rawDY_h[gammaNPos][ptPos]->Fill(MDzeroUPC->Dy->at(dI)); + rawDY_h[gammaNPos][nDPtBins]->Fill(MDzeroUPC->Dy->at(dI)); } if(ptPos >= 0 && yPos >= 0){ - rawDM_h[ptPos][yPos]->Fill(MDzeroUPC->Dmass->at(dI)); - rawDM_h[nDPtBins][yPos]->Fill(MDzeroUPC->Dmass->at(dI)); - rawDM_h[ptPos][nDYBins]->Fill(MDzeroUPC->Dmass->at(dI)); - rawDM_h[nDPtBins][nDYBins]->Fill(MDzeroUPC->Dmass->at(dI)); + rawDM_h[gammaNPos][ptPos][yPos]->Fill(MDzeroUPC->Dmass->at(dI)); + rawDM_h[gammaNPos][nDPtBins][yPos]->Fill(MDzeroUPC->Dmass->at(dI)); + rawDM_h[gammaNPos][ptPos][nDYBins]->Fill(MDzeroUPC->Dmass->at(dI)); + rawDM_h[gammaNPos][nDPtBins][nDYBins]->Fill(MDzeroUPC->Dmass->at(dI)); } } } delete MDzeroUPC; - + inFile_p->Close(); delete inFile_p; + + + if(nCurrEntries > totalNEntries) break; }//complete file loop //outfile handling outFile_p->cd(); + cutPass_h->Write("", TObject::kOverwrite); + delete cutPass_h; - for(unsigned int yI = 0; yI < nDYBins+1; ++yI){ - rawDPt_h[yI]->Write("", TObject::kOverwrite); - delete rawDPt_h[yI]; - } + for(auto const& gammaNVal : nGammaBinsStrMap){ + for(unsigned int yI = 0; yI < nDYBins+1; ++yI){ + rawDPt_h[gammaNVal.second][yI]->Write("", TObject::kOverwrite); + delete rawDPt_h[gammaNVal.second][yI]; + } - for(unsigned int pI = 0; pI < nDPtBins+1; ++pI){ - rawDY_h[pI]->Write("", TObject::kOverwrite); - delete rawDY_h[pI]; - } + for(unsigned int pI = 0; pI < nDPtBins+1; ++pI){ + rawDY_h[gammaNVal.second][pI]->Write("", TObject::kOverwrite); + delete rawDY_h[gammaNVal.second][pI]; + } - for(Int_t pI = 0; pI < nDPtBins+1; ++pI){ - for(unsigned int yI = 0; yI < nDYBins+1; ++yI){ - rawDM_h[pI][yI]->Write("", TObject::kOverwrite); - delete rawDM_h[pI][yI]; + for(Int_t pI = 0; pI < nDPtBins+1; ++pI){ + for(unsigned int yI = 0; yI < nDYBins+1; ++yI){ + rawDM_h[gammaNVal.second][pI][yI]->Write("", TObject::kOverwrite); + delete rawDM_h[gammaNVal.second][pI][yI]; + } } } - + //Write out the config so there is a record of production config_p->Write("config", TObject::kOverwrite); delete config_p; - + outFile_p->Close(); delete outFile_p; - + std::cout << "DataMCComp complete! return 0" << std::endl; return 0; } diff --git a/MainAnalysis/20241106_DzeroUPCDataMCComp/configs/dataMCCompare.config b/MainAnalysis/20241106_DzeroUPCDataMCComp/configs/dataMCCompare.config index 120e7065..aed799be 100644 --- a/MainAnalysis/20241106_DzeroUPCDataMCComp/configs/dataMCCompare.config +++ b/MainAnalysis/20241106_DzeroUPCDataMCComp/configs/dataMCCompare.config @@ -1,10 +1,11 @@ SKIMFILENAME: /data/NewSkims23_24/20241106_SkimOldReco23sample_DataAll.root -DPTBINS: 0,2,5,8,12,20 +DPTBINS: 2,5,8,12,20 DYBINS: -2,-1,0,1,2 NDMBINS: 50 DMBINSLOW: 1.6 DMBINSHIGH: 2.1 -NENTRIESTOPRINT: 1000000 \ No newline at end of file +NENTRIESTOPRINT: 1000000 +NENTRIESTOPROCESS: 10000000 \ No newline at end of file