diff --git a/MainAnalysis/20241106_DzeroUPCDataMCComp/DataMCComp.cpp b/MainAnalysis/20241106_DzeroUPCDataMCComp/DataMCComp.cpp new file mode 100644 index 00000000..6e72e7a2 --- /dev/null +++ b/MainAnalysis/20241106_DzeroUPCDataMCComp/DataMCComp.cpp @@ -0,0 +1,350 @@ +//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" + +//Local CommonCode +#include "Messenger.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 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; +} + +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 - + 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 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"); + 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; + + //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"); + DzeroUPCTreeMessenger *MDzeroUPC = new DzeroUPCTreeMessenger(inFile_p, "Tree"); + totalNEntries += MDzeroUPC->GetEntries(); + 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 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"); + //Use UPC DZeroTree + DzeroUPCTreeMessenger *MDzeroUPC = new DzeroUPCTreeMessenger(inFile_p, "Tree"); + + //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; + + 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[gammaNPos][yPos]->Fill(MDzeroUPC->Dpt->at(dI)); + rawDPt_h[gammaNPos][nDYBins]->Fill(MDzeroUPC->Dpt->at(dI)); + } + + if(ptPos >= 0){ + 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[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(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[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[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; +} + +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..aed799be --- /dev/null +++ b/MainAnalysis/20241106_DzeroUPCDataMCComp/configs/dataMCCompare.config @@ -0,0 +1,11 @@ +SKIMFILENAME: /data/NewSkims23_24/20241106_SkimOldReco23sample_DataAll.root + +DPTBINS: 2,5,8,12,20 +DYBINS: -2,-1,0,1,2 + +NDMBINS: 50 +DMBINSLOW: 1.6 +DMBINSHIGH: 2.1 + +NENTRIESTOPRINT: 1000000 +NENTRIESTOPROCESS: 10000000 \ 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