Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2024-06-26 07:06:07 

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 Christopher Dilks, Connor Pecar, Duane Byer, Sanghwa Park, Matthew McEneaney, Brian Page
0003 
0004 #include "Analysis.h"
0005 
0006 ClassImp(Analysis)
0007 
0008 using std::cout;
0009 using std::cerr;
0010 using std::endl;
0011 
0012 // constructor
0013 Analysis::Analysis(
0014   TString configFileName_,
0015   TString outfilePrefix_
0016 )
0017   : configFileName(configFileName_)
0018   , outfilePrefix(outfilePrefix_)
0019   , reconMethod("")
0020   , finalStateID("")
0021   , eleBeamEn(10.0)
0022   , ionBeamEn(100.0)
0023   , crossingAngle(-25.0)
0024   , totalCrossSection(1e6)
0025 {
0026   // available variables for binning
0027   // - availableBinSchemes is a map from variable name to variable title
0028   // - try to avoid using underscores in the variable name (they are okay in the title);
0029   //   convention is camel case, starting lowercase 
0030   /* DIS */
0031   availableBinSchemes.insert({ "x",     "x"           });
0032   availableBinSchemes.insert({ "q2",    "Q^{2}"       });
0033   availableBinSchemes.insert({ "w",     "W"           });
0034   availableBinSchemes.insert({ "y",     "y"           });
0035   /* single hadron */
0036   availableBinSchemes.insert({ "p",     "p"           });
0037   availableBinSchemes.insert({ "eta",   "#eta"        });
0038   availableBinSchemes.insert({ "pt",    "p_{T}"       }); // transverse to q, in ion rest frame
0039   availableBinSchemes.insert({ "ptLab", "p_{T}^{lab}" }); // transverse to xy plane, in lab frame
0040   availableBinSchemes.insert({ "z",     "z"           });
0041   availableBinSchemes.insert({ "qT",    "q_{T}"       });
0042   availableBinSchemes.insert({ "qTq",   "q_{T}/Q"     });
0043   availableBinSchemes.insert({ "mX",    "M_{X}"       });
0044   availableBinSchemes.insert({ "xF",    "x_{F}"       });
0045   availableBinSchemes.insert({ "phiH",  "#phi_{h}"    });
0046   availableBinSchemes.insert({ "phiS",  "#phi_{S}"    });
0047   availableBinSchemes.insert({ "tSpin", "spin"        });
0048   availableBinSchemes.insert({ "lSpin", "spinL"       });
0049   /* jets */
0050 #ifndef EXCLUDE_DELPHES
0051   availableBinSchemes.insert({ "JetPT", "jet p_{T}" });
0052   availableBinSchemes.insert({ "JetZ",  "jet z"     });
0053   availableBinSchemes.insert({ "JetEta", "jet eta" });
0054   availableBinSchemes.insert({ "JetE", "jet energy" });
0055 #endif
0056 
0057   // available final states
0058   // - specify which final states you want to include using `AddFinalState(TString name)`
0059   // - if you specify none, default final state(s) will be chosen for you
0060   availableBinSchemes.insert({"finalState","finalState"});
0061   AddBinScheme("finalState");
0062   // - finalState name (ID) -> title
0063   finalStateToTitle.insert({ "pipTrack", "#pi^{+} tracks" });
0064   finalStateToTitle.insert({ "pimTrack", "#pi^{-} tracks" });
0065   finalStateToTitle.insert({ "KpTrack",  "K^{+} tracks"   });
0066   finalStateToTitle.insert({ "KmTrack",  "K^{-} tracks"   });
0067   finalStateToTitle.insert({ "pTrack",   "p^{+} tracks"   });
0068   // - PID -> finalState ID
0069   PIDtoFinalState.insert({  211,  "pipTrack" });
0070   PIDtoFinalState.insert({ -211,  "pimTrack" });
0071   PIDtoFinalState.insert({  321,  "KpTrack"  });
0072   PIDtoFinalState.insert({ -321,  "KmTrack"  });
0073   PIDtoFinalState.insert({  2212, "pTrack"   });
0074 
0075   // kinematics reconstruction methods
0076   // - choose one of these methods using `SetReconMethod(TString name)`
0077   // - if you specify none, a default method will be chosen
0078   reconMethodToTitle.insert({ "Ele",    "Electron method"        });
0079   reconMethodToTitle.insert({ "DA",     "Double Angle method"    });
0080   reconMethodToTitle.insert({ "JB",     "Jacquet-Blondel method" });
0081   reconMethodToTitle.insert({ "Mixed",  "Mixed method"           });
0082   reconMethodToTitle.insert({ "Sigma",  "Sigma method"           });
0083   reconMethodToTitle.insert({ "eSigma", "eSigma method"          });
0084 
0085   // output sets to include
0086   // - use these to turn on/off certain sets of variables
0087   // - the default settings are set here; override them at the macro level
0088   includeOutputSet.insert({ "inclusive",      true  }); // inclusive kinematics
0089   includeOutputSet.insert({ "1h",             true  }); // single hadron kinematics
0090   includeOutputSet.insert({ "jets",           false }); // jet kinematics
0091   includeOutputSet.insert({ "depolarization", false }); // depolarization factors & ratios
0092 
0093   // common settings defaults
0094   // - these settings can be set at the macro level
0095   verbose           = false;
0096   writeSidisTree   = false;
0097   writeHFSTree      = false;
0098   writeParticleTree = false;
0099 
0100   maxEvents    = 0;
0101   useBreitJets = false;
0102   errorCntMax  = 1000;
0103   jetAlg       = 0; // Default to kT Algorithm
0104   jetRad       = 0.8;
0105   jetMin       = 1.0; // Minimum Jet pT
0106   jetMatchDR   = 0.5; // Delta R between true and reco jet to be considered matched
0107 
0108   weightInclusive = std::make_unique<WeightsUniform>();
0109   weightTrack     = std::make_unique<WeightsUniform>();
0110   weightJet       = std::make_unique<WeightsUniform>();
0111 
0112   // miscellaneous
0113   infiles.clear();
0114   entriesTot = 0;
0115   errorCnt = 0;
0116 };
0117 
0118 
0119 // input files
0120 //------------------------------------
0121 // add a group of files
0122 void Analysis::AddFileGroup(
0123     std::vector<std::string> fileNames,
0124     Long64_t totalEntries,
0125     Double_t xs,
0126     Double_t Q2min,
0127     Double_t Q2max,
0128     Double_t manualWeight
0129     )
0130 {
0131   // print
0132   fmt::print("-> AddFileGroup:  crossSection={},  Q2range = {} to {}\n",
0133       xs,
0134       Q2min,
0135       Q2max>Q2min ? std::to_string(Q2max) : "inf"
0136       );
0137   for (auto fileName : fileNames) fmt::print("   - {}\n", fileName);
0138 
0139   // fill vectors
0140   infiles.push_back(fileNames);
0141   Q2xsecs.push_back(xs);
0142   Q2mins.push_back(Q2min);
0143   Q2maxs.push_back(Q2max);
0144   Q2entries.push_back(totalEntries);
0145   if (manualWeight!=0)
0146       Q2weights.push_back(manualWeight);
0147     
0148   // check if the cross section for each group decreases; this is preferred
0149   // to make sure that `GetEventQ2Idx` behaves correctly
0150   for (std::size_t idx = 0; idx+1 < Q2xsecs.size(); ++idx) {
0151     if (Q2xsecs[idx] < Q2xsecs[idx+1] ) {
0152       cerr << "WARNING: Cross-sections should strictly decrease with stricter Q2min cuts; re-arrange your config file" << endl;
0153       PrintStdVector(Q2xsecs,"   cross sections");
0154     }
0155   }
0156 
0157   // lookup table, for tree number -> vector indices
0158   inLookup.clear();
0159   std::size_t total = 0;
0160   for (std::size_t idx = 0; idx < infiles.size(); ++idx) {
0161     total += infiles[idx].size();
0162     inLookup.resize(total, idx);
0163   }
0164 }
0165 
0166 
0167 // prepare for the analysis
0168 //------------------------------------
0169 void Analysis::Prepare() {
0170 
0171   // parse config file --------------------
0172   std::ifstream fin(configFileName);
0173   std::string line;
0174   bool debugParser = false;
0175   fmt::print("[+++] PARSING CONFIG FILE {}\n",configFileName);
0176 
0177   // vars
0178   Double_t xsec;
0179   Double_t Q2min;
0180   Double_t Q2max;
0181   Double_t manualWeight;
0182   Long64_t numEvents;
0183   std::vector<std::string> fileNames;
0184   bool readingListOfFiles;
0185   bool endGroupCalled;
0186 
0187   // reset vars
0188   auto ResetVars = [&] () {
0189     xsec      = totalCrossSection;
0190     Q2min     = 1.0;
0191     Q2max     = 0.0;
0192     manualWeight    = 0.0;
0193     numEvents = -1;
0194     fileNames.clear();
0195     readingListOfFiles = false;
0196     endGroupCalled     = false;
0197   };
0198   ResetVars();
0199 
0200   /* lambda to add a list of files to this analysis; wraps `Analysis::AddFileGroup`,
0201    * and checks the files beforehand
0202    */
0203   auto AddFiles = [this, &fileNames, &xsec, &Q2min, &Q2max, &numEvents, &manualWeight] () {
0204     Long64_t entries = 0;
0205     if(numEvents<0) {
0206       // get the number of entries, and check the file
0207       for(auto fileName : fileNames) {
0208         auto file = TFile::Open(fileName.c_str());
0209         if (file==nullptr || file->IsZombie()) {
0210           fmt::print(stderr,"ERROR: Couldn't open input file '{}'\n",fileName);
0211           return;
0212         }
0213         TTree *tree = file->Get<TTree>("Delphes");                  // fastsim
0214         if (tree == nullptr) tree = file->Get<TTree>("events");     // ATHENA, ePIC
0215         if (tree == nullptr) tree = file->Get<TTree>("event_tree"); // ECCE
0216         if (tree == nullptr) {
0217           fmt::print(stderr,"ERROR: Couldn't find tree in file '{}'\n",fileName);
0218         }
0219         else entries += tree->GetEntries();
0220         file->Close();
0221         delete file;
0222       }
0223     }
0224     else entries = numEvents;
0225     // add the file group
0226     AddFileGroup(fileNames, entries, xsec, Q2min, Q2max, manualWeight);
0227   };
0228 
0229   // loop over lines
0230   while (std::getline(fin, line)) {
0231 
0232     // chomp comments
0233     auto lineChomped = line.substr(0,line.find("#"));
0234     if(debugParser) fmt::print("\nlineChomped = \"{}\"\n",lineChomped);
0235 
0236     // each line will be parsed as one of the following:
0237     enum parseAs_enum { kNone, kSetting, kRootFile };
0238     int parseAs = kNone;
0239 
0240     // tokenize the line
0241     std::string token, bufferKey, bufferVal;
0242     std::stringstream lineStream(lineChomped);
0243     while (std::getline(lineStream, token, ' ')) {
0244       // classify `parseAs` and fill buffers
0245       if (token.rfind(":",0)==0) { // parse as setting name
0246         parseAs   = kSetting;
0247         bufferKey = token;
0248       }
0249       else if (parseAs==kSetting) { // parse as setting value
0250         if (token!="" && bufferVal=="") bufferVal = token;
0251       }
0252       else { // parse as root file name
0253         parseAs   = kRootFile;
0254         bufferKey = token;
0255       }
0256       // if(debugParser) fmt::print("  token = \"{}\"\n",token);
0257     } // tokenize
0258 
0259     // parse buffers
0260     if (parseAs==kSetting) {
0261       // if we were reading a list of files, we're done with this list; add the files and reset
0262       if(readingListOfFiles || bufferKey==":endGroup") {
0263         if(debugParser) fmt::print("-> new setting, add current list of `fileNames`, then reset\n");
0264         AddFiles();
0265         ResetVars();
0266       }
0267       // parse setting value(s)
0268       if(debugParser) fmt::print("  parse as setting: bufferKey='{}' bufferVal='{}'\n",bufferKey,bufferVal);
0269       if(bufferVal=="" && bufferKey!=":endGroup")
0270         fmt::print(stderr,"ERROR: setting '{}' has no associated value\n",bufferKey);
0271       else if (bufferKey==":eleBeamEn")         eleBeamEn         = std::stod(bufferVal);
0272       else if (bufferKey==":ionBeamEn")         ionBeamEn         = std::stod(bufferVal);
0273       else if (bufferKey==":crossingAngle")     crossingAngle     = std::stod(bufferVal);
0274       else if (bufferKey==":totalCrossSection") totalCrossSection = std::stod(bufferVal);
0275       else if (bufferKey==":q2min")             Q2min             = std::stod(bufferVal);
0276       else if (bufferKey==":q2max")             Q2max             = std::stod(bufferVal);
0277       else if (bufferKey==":crossSection")      xsec              = std::stod(bufferVal);
0278       else if (bufferKey==":Weight")            manualWeight      = std::stod(bufferVal);
0279       else if (bufferKey==":numEvents")         numEvents         = std::stoll(bufferVal);
0280       else if (bufferKey==":endGroup")          endGroupCalled    = true;
0281       else
0282         fmt::print(stderr,"ERROR: unknown setting \"{}\"\n",bufferKey);
0283       if(debugParser) fmt::print("  setting:  \"{}\" = \"{}\"\n",bufferKey,bufferVal);
0284     }
0285     else if (parseAs==kRootFile) {
0286       readingListOfFiles = true;
0287       fileNames.push_back(bufferKey);
0288       if(debugParser) fmt::print("  rootFile: \"{}\"\n",bufferKey);
0289     }
0290 
0291   } // loop over lines
0292 
0293   // add last list of files
0294   if(debugParser) fmt::print("-> EOF; add last list of `fileNames`\n");
0295   if(!endGroupCalled) AddFiles();
0296   fmt::print("[+++] PARSING COMPLETE\n\n");
0297 
0298   if (infiles.empty()) {
0299     cerr << "ERROR: no input files have been specified" << endl;
0300     return;
0301   }
0302 
0303   // print configuration ----------------------------------
0304   fmt::print("{:=<50}\n","CONFIGURATION: ");
0305   fmt::print("{:>30} = {} GeV\n",  "eleBeamEn",         eleBeamEn);
0306   fmt::print("{:>30} = {} GeV\n",  "ionBeamEn",         ionBeamEn);
0307   fmt::print("{:>30} = {} mrad\n", "crossingAngle",     crossingAngle);
0308   fmt::print("{:>30} = {}\n",      "totalCrossSection", totalCrossSection);
0309   fmt::print("{:>30} = {}\n",      "reconMethod",       reconMethod);
0310   fmt::print("{:-<50}\n","");
0311   PrintStdVector(Q2mins,"Q2mins");
0312   PrintStdVector(Q2maxs,"Q2maxs");
0313   PrintStdVector(Q2xsecs,"Q2xsecs");
0314   PrintStdVector(Q2entries,"Q2entries");
0315   // PrintStdVector(inLookup,"inLookup");
0316   // fmt::print("{:-<50}\n",""); PrintStdMap(availableBinSchemes,"availableBinSchemes");
0317   fmt::print("{:-<50}\n",""); PrintStdMap(includeOutputSet,"includeOutputSet");
0318   fmt::print("{:=<50}\n","");
0319 
0320   // set output file name
0321   outfileName = "out/"+outfilePrefix+".root";
0322 
0323   // open output file
0324   cout << "-- output file: " << outfileName << endl;
0325   outFile = new TFile(outfileName,"RECREATE");
0326 
0327   // instantiate shared objects
0328   kin        = std::make_shared<Kinematics>(eleBeamEn,ionBeamEn,crossingAngle);
0329   kinTrue    = std::make_shared<Kinematics>(eleBeamEn, ionBeamEn, crossingAngle);
0330 #ifndef EXCLUDE_DELPHES
0331   kinJet     = std::make_shared<KinematicsJets>(eleBeamEn, ionBeamEn, crossingAngle);
0332   kinJetTrue = std::make_shared<KinematicsJets>(eleBeamEn, ionBeamEn, crossingAngle);
0333 #endif
0334   ST         = std::make_unique<SidisTree>("tree",kin,kinTrue);
0335   HFST       = std::make_unique<HFSTree>("hfstree",kin,kinTrue);
0336   PT         = std::make_unique<ParticleTree>("ptree");
0337 
0338   // if including jets, define a `jet` final state
0339 #ifndef EXCLUDE_DELPHES
0340   if(includeOutputSet["jets"]) {
0341     finalStateToTitle.insert({ "jet", "jets" });
0342     AddFinalState("jet");
0343   }
0344 #endif
0345 
0346   // determine if only including `inclusive` output set
0347   includeOutputSet.insert({ "inclusive_only",
0348       includeOutputSet["inclusive"]
0349       && !includeOutputSet["1h"]
0350       && !includeOutputSet["jets"]
0351       });
0352 
0353   // if there are no final states defined, default to definitions here:
0354   if(BinScheme("finalState")->GetNumBins()==0) {
0355     std::cout << "NOTE: adding pi+ tracks for final state, since you specified none" << std::endl;
0356     AddFinalState("pipTrack");
0357   };
0358 
0359   // if no reconstruction method is set, choose a default here
0360   if(reconMethod=="") {
0361     std::cout << "NOTE: no recon method specified, default to electron method" << std::endl;
0362     SetReconMethod("Ele");
0363   };
0364 
0365   // build HistosDAG with specified binning
0366   HD = std::make_shared<HistosDAG>();
0367   HD->Build(binSchemes);
0368 
0369   // tell HD what values to associate for each BinScheme name
0370   // - these are the values that will be used for CutDef checks
0371   // - the lambda must return a Double_t
0372   /* DIS */
0373   HD->SetBinSchemeValue("x",     [this](){ return kin->x;                       });
0374   HD->SetBinSchemeValue("q2",    [this](){ return kin->Q2;                      });
0375   HD->SetBinSchemeValue("w",     [this](){ return kin->W;                       });
0376   HD->SetBinSchemeValue("y",     [this](){ return kin->y;                       });
0377   /* single hadron */
0378   HD->SetBinSchemeValue("p",     [this](){ return kin->pLab;                    });
0379   HD->SetBinSchemeValue("eta",   [this](){ return kin->etaLab;                  });
0380   HD->SetBinSchemeValue("pt",    [this](){ return kin->pT;                      });
0381   HD->SetBinSchemeValue("ptLab", [this](){ return kin->pTlab;                   });
0382   HD->SetBinSchemeValue("z",     [this](){ return kin->z;                       });
0383   HD->SetBinSchemeValue("qT",    [this](){ return kin->qT;                      });
0384   HD->SetBinSchemeValue("qTq",   [this](){ return kin->qT/TMath::Sqrt(kin->Q2); });
0385   HD->SetBinSchemeValue("mX",    [this](){ return kin->mX;                      });
0386   HD->SetBinSchemeValue("xF",    [this](){ return kin->xF;                      });
0387   HD->SetBinSchemeValue("phiH",  [this](){ return kin->phiH;                    });
0388   HD->SetBinSchemeValue("phiS",  [this](){ return kin->phiS;                    });
0389   HD->SetBinSchemeValue("tSpin", [this](){ return (Double_t)kin->tSpin;         });
0390   HD->SetBinSchemeValue("lSpin", [this](){ return (Double_t)kin->lSpin;         });
0391   /* jets */
0392 #ifndef EXCLUDE_DELPHES
0393   HD->SetBinSchemeValue("JetPT", [this](){ return kinJet->pTjet;                });
0394   HD->SetBinSchemeValue("JetZ",  [this](){ return kinJet->zjet;                 });
0395   HD->SetBinSchemeValue("JetEta", [this](){ return kinJet->etajet;              });
0396   HD->SetBinSchemeValue("JetE", [this](){ return kinJet->ejet;                  });
0397 #endif
0398 
0399   // some bin schemes values are checked here, instead of by CutDef checks ("External" cut type)
0400   // - the lambda must return a boolean
0401   HD->SetBinSchemeValueExternal("finalState", [this](Node *N){ return N->GetCut()->GetCutID() == finalStateID; });
0402 
0403 
0404   // DEFINE HISTOGRAMS ------------------------------------
0405   // - whether they are defined is controlled by `includeOutputSet` settings
0406   // - `Histos::Hist` calls should check for existence
0407   HD->Payload([this](Histos *HS){
0408     // -- inclusive kinematics
0409     if(includeOutputSet["inclusive"]) {
0410       // -- Full phase space histogram
0411       HS->DefineHist4D(
0412           "full_xsec",
0413           "x","Q^{2}","z","p_{T}",
0414           "","GeV^{2}","","GeV",
0415           NBINS_FULL,1e-3,1,
0416           NBINS_FULL,1,100,
0417           NBINS_FULL,0,1,
0418           NBINS_FULL,0,2,
0419           true,true
0420           );
0421       // -- DIS kinematics
0422       HS->DefineHist2D("Q2vsX","x","Q^{2}","","GeV^{2}",
0423           NBINS,1e-3,1,
0424           NBINS,1,3000,
0425           true,true
0426           );
0427       HS->DefineHist1D("Q2","Q2","GeV",NBINS,1.0,3000,true,true);
0428       HS->DefineHist1D("x","x","",NBINS,1e-3,1.0,true,true);
0429       HS->DefineHist1D("y","y","",NBINS,1e-3,1,true);
0430       HS->DefineHist1D("W","W","GeV",NBINS,0,50);
0431     }
0432 
0433     // -- single hadron kinematics
0434     if(includeOutputSet["1h"]) {
0435       // -- hadron 4-momentum
0436       HS->DefineHist1D("pLab","p_{lab}","GeV",NBINS,0,10);
0437       HS->DefineHist1D("pTlab","p_{T}^{lab}","GeV",NBINS,1e-2,3,true);
0438       HS->DefineHist1D("etaLab","#eta_{lab}","",NBINS,-5,5);
0439       HS->DefineHist1D("phiLab","#phi_{lab}","",NBINS,-TMath::Pi(),TMath::Pi());
0440       // -- hadron kinematics
0441       HS->DefineHist1D("z","z","",NBINS,0,1);
0442       HS->DefineHist1D("pT","p_{T}","GeV",NBINS,1e-2,3,true);
0443       HS->DefineHist1D("qT","q_{T}","GeV",NBINS,1e-2,5,true);
0444       HS->DefineHist1D("qTq","q_{T}/Q","",NBINS,1e-2,3,true);
0445       HS->DefineHist1D("mX","m_{X}","GeV",NBINS,0,40);
0446       HS->DefineHist1D("phiH","#phi_{h}","",NBINS,-TMath::Pi(),TMath::Pi());
0447       HS->DefineHist1D("phiS","#phi_{S}","",NBINS,-TMath::Pi(),TMath::Pi());
0448       HS->DefineHist2D("phiHvsPhiS","#phi_{S}","#phi_{h}","","",
0449           25,-TMath::Pi(),TMath::Pi(),
0450           25,-TMath::Pi(),TMath::Pi());
0451       HS->DefineHist1D("phiSivers","#phi_{Sivers}","",NBINS,-TMath::Pi(),TMath::Pi());
0452       HS->DefineHist1D("phiCollins","#phi_{Collins}","",NBINS,-TMath::Pi(),TMath::Pi());
0453       HS->DefineHist2D("etaVsP","p","#eta","GeV","",
0454           NBINS,0.1,100,
0455           NBINS,-4,4,
0456           true,false,true
0457           );
0458       Double_t etabinsCoarse[] = {-4.0,-1.0,1.0,4.0};
0459       Double_t pbinsCoarse[] = {0.1,1,10,100};
0460       HS->DefineHist2D("etaVsPcoarse","p","#eta","GeV","",
0461           3, pbinsCoarse,
0462           3, etabinsCoarse,
0463           true,false
0464           );
0465       // -- single-hadron cross sections
0466       //HS->DefineHist1D("Q_xsec","Q","GeV",10,0.5,10.5,false,true); // linear
0467       HS->DefineHist1D("Q_xsec","Q","GeV",NBINS,1.0,3000,true,true); // log
0468       HS->Hist("Q_xsec")->SetMinimum(1e-10);
0469       // -- resolutions
0470       HS->DefineHist1D("x_Res",  "x-x_{true}/x_{true}",             "", NBINS, -1.0, 1.0);
0471       HS->DefineHist1D("y_Res",  "y-y_{true}/y_{true}",             "", NBINS, -1.0, 1.0);
0472       HS->DefineHist1D("Q2_Res", "Q^{2}-Q^{2}_{true}/Q^{2}_{true}", "", NBINS, -1.0, 1.0);
0473       HS->DefineHist1D("W_Res",  "W-W_{true}/W_{true}",             "", NBINS, -1.0, 1.0);
0474       HS->DefineHist1D("Nu_Res", "#nu-#nu_{true}/#nu_{true}",       "", NBINS, -1.0, 1.0);
0475       HS->DefineHist1D("pT_Res", "p_{T}-p_{T}^{true}/p_{T}^{true}", "", NBINS, -1.0, 1.0);
0476       HS->DefineHist1D("z_Res",  "z-z_{true}/z_{true}",             "", NBINS, -1.0, 1.0);
0477       HS->DefineHist1D("mX_Res", "m_{X}-m_{X}^{true}/m_{X}^{true}", "", NBINS, -1.0, 1.0);
0478       HS->DefineHist1D("xF_Res", "x_{F}-x_{F}^{true}/x_{F}^{true}", "", NBINS, -1.0, 1.0);
0479       HS->DefineHist1D("phiH_Res", "#phi_{h}-#phi_{h}^{true}", "", NBINS, -TMath::Pi(), TMath::Pi()); // absolute resolution
0480       HS->DefineHist1D("phiS_Res", "#phi_{S}-#phi_{S}^{true}", "", NBINS, -TMath::Pi(), TMath::Pi()); // absolute resolution
0481       HS->DefineHist2D("Q2vsXtrue","x","Q^{2}","","GeV^{2}",
0482           20,1e-4,1,
0483           10,1,1e4,
0484           true,true
0485           );
0486       HS->DefineHist2D("Q2vsXpurity","x","Q^{2}","","GeV^{2}",
0487           20,1e-4,1,
0488           10,1,1e4,
0489           true,true
0490           );
0491       HS->DefineHist2D("Q2vsX_zres","x","Q^{2}","","GeV^{2}",
0492           20,1e-4,1,
0493           10,1,1e4,
0494           true,true
0495           );
0496       HS->DefineHist2D("Q2vsX_pTres","x","Q^{2}","","GeV^{2}",
0497           20,1e-4,1,
0498           10,1,1e4,
0499           true,true
0500           );
0501       HS->DefineHist2D("Q2vsX_phiHres","x","Q^{2}","","GeV^{2}",
0502           20,1e-4,1,
0503           10,1,1e4,
0504           true,true
0505           );
0506       // -- reconstructed vs. generated
0507       HS->DefineHist2D("x_RvG","generated x","reconstructed x","","",
0508           NBINS,1e-3,1,
0509           NBINS,1e-3,1,
0510           true,true
0511           );
0512       HS->DefineHist2D("phiH_RvG","generated #phi_{h}","reconstructed #phi_{h}","","",
0513           NBINS,-TMath::Pi(),TMath::Pi(),
0514           NBINS,-TMath::Pi(),TMath::Pi()
0515           );
0516       HS->DefineHist2D("phiS_RvG","generated #phi_{S}","reconstructed #phi_{S}","","",
0517           NBINS,-TMath::Pi(),TMath::Pi(),
0518           NBINS,-TMath::Pi(),TMath::Pi()
0519           );
0520     }
0521 
0522     // -- jet kinematics
0523 #ifndef EXCLUDE_DELPHES
0524     if(includeOutputSet["jets"]) {
0525       HS->DefineHist1D("JetPT",    "jet p_{T}",      "GeV", NBINS, 1e-2, 50,  true, true);
0526       HS->DefineHist1D("JetMT",    "jet m_{T}",      "GeV", NBINS, 1e-2, 20,  true, true);
0527       HS->DefineHist1D("JetZ",     "jet z",          "",    NBINS, 0,    1);
0528       HS->DefineHist1D("JetEta",   "jet #eta_{lab}", "",    NBINS, -5,   5);
0529       HS->DefineHist1D("JetQT",    "jet q_{T}",      "GeV", NBINS, 0,    10.0);
0530       HS->DefineHist1D("JetJperp", "j_{#perp}",      "GeV", NBINS, 1e-2, 3.0, true, true);
0531       HS->DefineHist1D("JetQTQ",   "jet q_{T}/Q",    "",    NBINS, 0,    3.0);
0532       HS->DefineHist1D("JetE","jet Energy","GeV", NBINS, 1e-2, 100);
0533       HS->DefineHist1D("JetM","jet Mass","GeV", NBINS, 1e-2, 20);
0534 
0535       HS->DefineHist2D("JetPTVsEta","Eta","pT","GeV","",100,-5,5,100,0,50,false,false);
0536 
0537       // Resolution Plos
0538       HS->DefineHist1D("JetDeltaR","deltaR between matched reco and truth jet","",1000,-2.,10.);
0539  
0540       HS->DefineHist2D("JetPTTrueVsReco","Reco pT","True pT","GeV","GeV",100,0.,50.,100,0.,50.,false,false);
0541       HS->DefineHist2D("JetETrueVsReco","Reco E","True E","GeV","GeV",100,0.,100.,100,0.,100.,false,false);
0542  
0543       HS->DefineHist2D("JetResPTVsTruePT","True pT","(Reco-True)/True pT","GeV","",100,0.,50.,10000,-10.,10.,false,false);
0544       HS->DefineHist2D("JetResEVsTrueE","True E","(Reco-True)/True E","GeV","",100,0.,100.,10000,-10.,10.,false,false);
0545       
0546       HS->DefineHist2D("JetResPTVsRecoPT","Reco pT","(Reco-True)/True pT","GeV","",100,0.,50.,10000,-10.,10.,false,false);
0547       HS->DefineHist2D("JetResEVsRecoE","Reco E","(Reco-True)/True E","GeV","",100,0.,100.,10000,-10.,10.,false,false);
0548     }
0549 #endif
0550 
0551     // -- depolarization
0552     if(includeOutputSet["depolarization"]) {
0553       std::map<TString,TString> depols = {
0554         { "epsilon", "#epsilon" },
0555         { "depolA",  "A"        },
0556         { "depolB",  "B"        },
0557         { "depolV",  "V"        },
0558         { "depolBA", "B/A"      },
0559         { "depolCA", "C/A"      },
0560         { "depolVA", "V/A"      },
0561         { "depolWA", "W/A"      }
0562       };
0563       for(auto [name,title] : depols) {
0564         HS->DefineHist2D(name+"vsQ2", "Q^{2}", title, "GeV^{2}", "", NBINS, 1,    3000, NBINS, 0, 2.5, true, false);
0565         HS->DefineHist2D(name+"vsY",  "y",     title, "",        "", NBINS, 5e-3, 1,    NBINS, 0, 2.5, true, false);
0566         HS->DefineHist3D(name+"vsQ2vsX",
0567             "x",   "Q^{2}",   title,
0568             "",    "GeV^{2}", "",
0569             NBINS, 1e-3,      1,
0570             NBINS, 1,         3000,
0571             NBINS, 0,         2.5,
0572             true,  true,      false
0573             );
0574       }
0575     }
0576 
0577   });
0578   HD->ExecuteAndClearOps();
0579 
0580 
0581   // initialize total weights
0582   wInclusiveTotal = 0.;
0583   wTrackTotal     = 0.;
0584   wJetTotal       = 0.;
0585 };
0586 
0587 void Analysis::CalculateEventQ2Weights() {
0588   Q2weights.resize(Q2xsecs.size());
0589   entriesTot = 0;
0590   for (Long64_t entry : Q2entries) {
0591     entriesTot += entry;
0592   }
0593   cout << "Q2 weighting info:" << endl;
0594   for (Int_t idx = 0; idx < Q2xsecs.size(); ++idx) {
0595     // calculate total luminosity, and the luminosity that contains this Q2 range
0596     Double_t lumiTotal = Double_t(entriesTot)     / totalCrossSection;
0597     Double_t lumiThis  = Double_t(Q2entries[idx]) / Q2xsecs[idx];
0598     //// alternative `lumiThis`: try to correct for overlapping Q2 ranges
0599     lumiThis = 0.; // reset
0600     for (Int_t j = 0; j < Q2xsecs.size(); ++j) {
0601       // check if Q2 range `j` contains the Q2 range `idx`; if so, include its luminosity
0602       if (InQ2Range(Q2mins[idx], Q2mins[j], Q2maxs[j]) &&
0603           InQ2Range(Q2maxs[idx], Q2mins[j], Q2maxs[j], true))
0604       {
0605         lumiThis += Double_t(Q2entries[j]) / Q2xsecs[j];
0606       }
0607     }
0608     // calculate the weight for this Q2 range
0609     if(Q2weights[idx]==0)
0610         Q2weights[idx] = lumiTotal / lumiThis;
0611     cout << "\tQ2 > "     << Q2mins[idx];
0612     if(Q2maxs[idx]>0) cout << " && Q2 < " << Q2maxs[idx];
0613     cout << ":" << endl;
0614     cout << "\t\tcount    = " << Q2entries[idx] << endl;
0615     cout << "\t\txsec     = " << Q2xsecs[idx] << endl;
0616     cout << "\t\tweight   = " << Q2weights[idx] << endl;
0617   }
0618 }
0619 
0620 Double_t Analysis::GetEventQ2Weight(Double_t Q2, Int_t guess) {
0621   Int_t idx = GetEventQ2Idx(Q2, guess);
0622   if (idx == -1) {
0623     return 0.;
0624   } else {
0625     return Q2weights[idx];
0626   }
0627 }
0628 
0629 // get the `idx` for this value of `Q2`; tries to find the most restrictive Q2 range
0630 // that contains the given `Q2` value, in order to assign the correct weight
0631 Int_t Analysis::GetEventQ2Idx(Double_t Q2, Int_t guess) {
0632   // return guess; // just use the Q2 range, according to which ROOT file the event came from
0633   Int_t idx = -1;
0634   for (int k=0; k<Q2mins.size(); k++) {
0635     if (InQ2Range(Q2,Q2mins[k],Q2maxs[k])) idx = k;
0636   }
0637   if (idx<0) {
0638     // fmt::print(stderr,"WARNING: Q2={} not in any Q2 range\n",Q2); // usually just below the smallest Q2 min
0639     idx = 0; // assume the least-restrictive range
0640   }
0641   else if(idx<guess){
0642       // When crossingAngle!=0, some event can be reconstructed with trueQ2 < Q2min of the Monte Carlo.
0643       // If so, set the idx=guess, i.e. just assume the event was generated in its file's Q2range
0644       idx = guess; 
0645   }
0646   return idx;
0647 }
0648 
0649 
0650 // finish the analysis
0651 //-----------------------------------
0652 void Analysis::Finish() {
0653 
0654   // reset HD, to clean up after the event loop
0655   HD->ActivateAllNodes();
0656   HD->ClearOps();
0657 
0658   // calculate integrated luminosity
0659   Double_t lumi = wTrackTotal/totalCrossSection; // [nb^-1]
0660   cout << "Integrated Luminosity:       " << lumi << "/nb" << endl;
0661   cout << sep << endl;
0662 
0663   // calculate cross sections, and print yields
0664   HD->Initial([this](){ cout << sep << endl << "Histogram Entries:" << endl; });
0665   HD->Final([this](){ cout << sep << endl; });
0666   HD->Payload([&lumi](Histos *H){
0667     auto h_Q2vsX       = H->Hist("Q2vsX",true);
0668     auto h_Q_xsec      = H->Hist("Q_xsec",true);
0669     auto h_Q2vsXpurity = H->Hist("Q2vsXpurity",true);
0670     if(h_Q2vsX) {
0671       cout << H->GetSetTitle() << " ::: "
0672            << h_Q2vsX->GetEntries()
0673            << endl;
0674     }
0675     // calculate cross sections
0676     if(h_Q_xsec) h_Q_xsec->Scale(1./lumi); // TODO: generalize (`if (name contains "xsec") ...`)
0677     // divide resolution plots by true counts per x-Q2 bin
0678     if(h_Q2vsXpurity) {
0679       H->Hist("Q2vsXpurity")->Divide(H->Hist("Q2vsXtrue"));
0680       H->Hist("Q2vsX_zres")->Divide(H->Hist("Q2vsXtrue"));
0681       H->Hist("Q2vsX_pTres")->Divide(H->Hist("Q2vsXtrue"));
0682       H->Hist("Q2vsX_phiHres")->Divide(H->Hist("Q2vsXtrue"));        
0683     }
0684   });
0685   HD->ExecuteAndClearOps();
0686 
0687   // write histograms
0688   cout << sep << endl;
0689   cout << "writing ROOT file..." << endl;
0690   outFile->cd();
0691   if(writeSidisTree)   ST->WriteTree();
0692   if(writeHFSTree)      HFST->WriteTree();
0693   if(writeParticleTree) PT->WriteTree();
0694   HD->Payload([this](Histos *H){ H->WriteHists(outFile); }); HD->ExecuteAndClearOps();
0695   HD->Payload([this](Histos *H){ H->Write(); }); HD->ExecuteAndClearOps();
0696   std::vector<Double_t> vec_wInclusiveTotal { wInclusiveTotal };
0697   std::vector<Double_t> vec_wTrackTotal     { wTrackTotal     };
0698   std::vector<Double_t> vec_wJetTotal       { wJetTotal       };
0699   outFile->WriteObject(&Q2xsecs, "XsTotal");
0700   outFile->WriteObject(&vec_wInclusiveTotal, "WeightInclusiveTotal");
0701   outFile->WriteObject(&vec_wTrackTotal,     "WeightTrackTotal");
0702   outFile->WriteObject(&vec_wJetTotal,       "WeightJetTotal");
0703   outFile->WriteObject(&total_events,    "TotalEvents");
0704     
0705   // write binning schemes
0706   for(auto const &kv : binSchemes) {
0707     if(kv.second->GetNumBins()>0) kv.second->Write("binset__"+kv.first);
0708   };
0709 
0710   // close output
0711   outFile->Close();
0712   cout << outfileName << " written." << endl;
0713 };
0714 
0715 
0716 // access bin scheme by name
0717 //------------------------------------
0718 BinSet *Analysis::BinScheme(TString varname) {
0719   BinSet *ret;
0720   try { ret = binSchemes.at(varname); }
0721   catch(const std::out_of_range &ex) {
0722     cerr << "ERROR: bin scheme "
0723          << varname << " not found" << endl;
0724     return nullptr;
0725   };
0726   return ret;
0727 };
0728 
0729 
0730 // add a new bin scheme
0731 //------------------------------------
0732 void Analysis::AddBinScheme(TString varname) {
0733   TString vartitle;
0734   // TODO [low priority]: would be nice to make this lookup case insensitive
0735   try { vartitle = availableBinSchemes.at(varname); }
0736   catch(const std::out_of_range &ex) {
0737     cerr << "ERROR: bin scheme "
0738          << varname << " not available... skipping..." << endl;
0739     return;
0740   };
0741   if(binSchemes.find(varname)==binSchemes.end()) { // (duplicate prevention)
0742     BinSet *B = new BinSet(varname,vartitle);
0743     binSchemes.insert({varname,B});
0744   };
0745 };
0746 
0747 
0748 // add a final state bin
0749 //------------------------------------
0750 void Analysis::AddFinalState(TString finalStateN) {
0751   TString finalStateT;
0752   try { finalStateT = finalStateToTitle.at(finalStateN); }
0753   catch(const std::out_of_range &ex) {
0754     cerr << "ERROR: final state "
0755          << finalStateN << " not available... skipping..." << endl;
0756     return;
0757   };
0758   BinScheme("finalState")->BuildExternalBin(finalStateN,finalStateT);
0759   activeFinalStates.insert(finalStateN);
0760   fmt::print("AddFinalState: name='{}'\n               title='{}'\n",finalStateN,finalStateT);
0761 };
0762 
0763 
0764 // FillHistos methods: check bins and fill associated histograms
0765 // - checks which bins the track/jet/etc. falls in
0766 // - fills the histograms in the associated Histos objects
0767 //--------------------------------------------------------------------------
0768 
0769 // fill histograms (called by other specific FillHistos* methods)
0770 void Analysis::FillHistos(std::function<void(Histos*)> fill_payload) {
0771   // check which bins to fill, and activate the ones for which all defined cuts pass
0772   // (activates their corresponding `HistosDAG` nodes)
0773   HD->CheckBins();
0774   // fill histograms, for activated bins only
0775   HD->Payload(fill_payload);
0776   // execute the payload
0777   // - save time and don't call `ClearOps` (next loop will overwrite lambda)
0778   // - called with `activeNodesOnly==true` since we only want to fill bins associated
0779   //   with this track
0780   HD->ExecuteOps(true);
0781 };
0782 
0783 // fill inclusive histograms
0784 void Analysis::FillHistosInclusive(Double_t wgt) {
0785   auto fill_payload = [this,wgt] (Histos *H) {
0786     H->FillHist2D("Q2vsX", kin->x,  kin->Q2, wgt);
0787     H->FillHist1D("Q2",    kin->Q2, wgt);
0788     H->FillHist1D("x",     kin->x,  wgt);
0789     H->FillHist1D("W",     kin->W,  wgt);
0790     H->FillHist1D("y",     kin->y,  wgt);
0791   };
0792   FillHistos(fill_payload);
0793 };
0794 
0795 // fill 1h (single-hadron) histograms
0796 void Analysis::FillHistos1h(Double_t wgt) {
0797   auto fill_payload = [this,wgt] (Histos *H) {
0798     // Full phase space.
0799     H->FillHist4D("full_xsec", kin->x, kin->Q2, kin->pT, kin->z, wgt);
0800     // hadron 4-momentum
0801     H->FillHist1D("pLab",   kin->pLab,   wgt);
0802     H->FillHist1D("pTlab",  kin->pTlab,  wgt);
0803     H->FillHist1D("etaLab", kin->etaLab, wgt);
0804     H->FillHist1D("phiLab", kin->phiLab, wgt);
0805     // hadron kinematics
0806     H->FillHist1D("z",  kin->z,  wgt);
0807     H->FillHist1D("pT", kin->pT, wgt);
0808     H->FillHist1D("qT", kin->qT, wgt);
0809     if(kin->Q2!=0) H->FillHist1D("qTq",kin->qT/TMath::Sqrt(kin->Q2),wgt);
0810     H->FillHist1D("mX",         kin->mX,   wgt);
0811     H->FillHist1D("phiH",       kin->phiH, wgt);
0812     H->FillHist1D("phiS",       kin->phiS, wgt);
0813     H->FillHist2D("phiHvsPhiS", kin->phiS, kin->phiH, wgt);
0814     H->FillHist1D("phiSivers",  Kinematics::AdjAngle(kin->phiH - kin->phiS), wgt);
0815     H->FillHist1D("phiCollins", Kinematics::AdjAngle(kin->phiH + kin->phiS), wgt);
0816     H->FillHist2D("etaVsP",       kin->pLab, kin->etaLab, wgt); // TODO: lab-frame p, or some other frame?
0817     H->FillHist2D("etaVsPcoarse", kin->pLab, kin->etaLab, wgt);
0818     // depolarization
0819     if(includeOutputSet["depolarization"]) { // not necessary, but done for performance
0820       std::map<TString,Double_t> depols = {
0821         { "epsilon", kin->epsilon },
0822         { "depolA",  kin->depolA  },
0823         { "depolB",  kin->depolB  },
0824         { "depolV",  kin->depolV  },
0825         { "depolBA", kin->depolP1 },
0826         { "depolCA", kin->depolP2 },
0827         { "depolVA", kin->depolP3 },
0828         { "depolWA", kin->depolP4 }
0829       };
0830       for(auto [name,val] : depols) {
0831         H->FillHist2D(name+"vsQ2",    kin->Q2, val,     wgt);
0832         H->FillHist2D(name+"vsY",     kin->y,  val,     wgt);
0833         H->FillHist3D(name+"vsQ2vsX", kin->x,  kin->Q2, val, wgt);
0834       }
0835     }
0836     // cross sections (divide by lumi after all events processed)
0837     H->FillHist1D("Q_xsec", TMath::Sqrt(kin->Q2), wgt);
0838     // resolutions
0839     auto fillRelativeResolution = [&H,&wgt] (auto name, auto rec, auto gen) { if(gen!=0) H->FillHist1D(name, (rec-gen)/gen, wgt); };
0840     fillRelativeResolution("x_Res",  kin->x,  kinTrue->x);
0841     fillRelativeResolution("y_Res",  kin->y,  kinTrue->y);
0842     fillRelativeResolution("Q2_Res", kin->Q2, kinTrue->Q2);
0843     fillRelativeResolution("W_Res",  kin->W,  kinTrue->W);
0844     fillRelativeResolution("Nu_Res", kin->Nu, kinTrue->Nu);
0845     fillRelativeResolution("pT_Res", kin->pT, kinTrue->pT);
0846     fillRelativeResolution("z_Res",  kin->z,  kinTrue->z);
0847     fillRelativeResolution("mX_Res", kin->mX, kinTrue->mX);
0848     fillRelativeResolution("xF_Res", kin->xF, kinTrue->xF);
0849     H->FillHist1D("phiH_Res",  Kinematics::AdjAngle(kin->phiH - kinTrue->phiH), wgt );
0850     H->FillHist1D("phiS_Res",  Kinematics::AdjAngle(kin->phiS - kinTrue->phiS), wgt );
0851     H->FillHist2D("Q2vsXtrue", kinTrue->x,            kinTrue->Q2, wgt);
0852     if(kinTrue->z!=0)
0853       H->FillHist2D("Q2vsX_zres", kinTrue->x, kinTrue->Q2, wgt*( fabs(kinTrue->z - kin->z)/(kinTrue->z) ) );
0854     if(kinTrue->pT!=0)
0855       H->FillHist2D("Q2vsX_pTres", kinTrue->x, kinTrue->Q2, wgt*( fabs(kinTrue->pT - kin->pT)/(kinTrue->pT) ) );
0856     H->FillHist2D("Q2vsX_phiHres", kinTrue->x, kinTrue->Q2, wgt*( fabs(Kinematics::AdjAngle(kinTrue->phiH - kin->phiH) ) ) );
0857 
0858     auto htrue = H->Hist("Q2vsXtrue",true);
0859     if(htrue!=nullptr) {
0860       if( htrue->FindBin(kinTrue->x,kinTrue->Q2) == htrue->FindBin(kin->x,kin->Q2) )
0861         H->FillHist2D("Q2vsXpurity", kin->x, kin->Q2, wgt);
0862     }
0863     // -- reconstructed vs. generated
0864     H->FillHist2D("x_RvG",    kinTrue->x,    kin->x,    wgt);
0865     H->FillHist2D("phiH_RvG", kinTrue->phiH, kin->phiH, wgt);
0866     H->FillHist2D("phiS_RvG", kinTrue->phiS, kin->phiS, wgt);
0867   };
0868   FillHistos(fill_payload);
0869 };
0870 
0871 // fill jet histograms
0872 void Analysis::FillHistosJets(Double_t wgt) {
0873 #ifndef EXCLUDE_DELPHES
0874   auto fill_payload = [this,wgt] (Histos *H) {
0875     // jet kinematics
0876     H->FillHist1D("JetPT",  kinJet->pTjet,  wgt);
0877     H->FillHist1D("JetMT",  kinJet->mTjet,  wgt);
0878     H->FillHist1D("JetZ",   kinJet->zjet,   wgt);
0879     H->FillHist1D("JetEta", kinJet->etajet, wgt);
0880     H->FillHist1D("JetQT",  kinJet->qTjet,  wgt);
0881     H->FillHist1D("JetE",   kinJet->ejet,   wgt);
0882     H->FillHist1D("JetM",   kinJet->mjet,   wgt);
0883 
0884     H->FillHist2D("JetPTVsEta", kinJet->etajet, kinJet->pTjet, wgt);
0885 
0886     // Fill Resolution Histos
0887     H->FillHist1D("JetDeltaR", kinJet->deltaRjet, wgt);
0888 
0889     H->FillHist2D("JetPTTrueVsReco", kinJet->pTjet, kinJet->pTmtjet, wgt);
0890     H->FillHist2D("JetETrueVsReco", kinJet->ejet, kinJet->emtjet, wgt);
0891 
0892     H->FillHist2D("JetResPTVsTruePT", kinJet->pTmtjet, (kinJet->pTjet - kinJet->pTmtjet)/(kinJet->pTmtjet), wgt);
0893     H->FillHist2D("JetResEVsTrueE", kinJet->emtjet, (kinJet->ejet - kinJet->emtjet)/(kinJet->emtjet), wgt);
0894 
0895     H->FillHist2D("JetResPTVsRecoPT", kinJet->pTjet, (kinJet->pTjet - kinJet->pTmtjet)/(kinJet->pTmtjet), wgt);
0896     H->FillHist2D("JetResEVsRecoE", kinJet->ejet, (kinJet->ejet - kinJet->emtjet)/(kinJet->emtjet), wgt);
0897 
0898     if(kinJet->Q2!=0) H->FillHist1D("JetQTQ", kinJet->qTjet/sqrt(kinJet->Q2), wgt);
0899     for(int j = 0; j < kinJet->jperp.size(); j++) {
0900       H->FillHist1D("JetJperp", kinJet->jperp[j], wgt);
0901     };
0902   };
0903   FillHistos(fill_payload);
0904 #endif
0905 };
0906 
0907 
0908 // print an error; if more than `errorCntMax` errors are printed, printing is suppressed
0909 void Analysis::ErrorPrint(std::string message) {
0910   errorCnt++;
0911   if(errorCnt <= errorCntMax) fmt::print(stderr,"{}\n",message);
0912   if(errorCnt == errorCntMax) fmt::print(stderr,"... {} errors printed; suppressing the rest ...\n",errorCnt);
0913 }
0914 
0915 // destructor
0916 Analysis::~Analysis() {
0917   for(auto it : binSchemes)
0918     if(it.second) delete it.second;
0919   binSchemes.clear();
0920 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.5 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!