diff --git a/macro/dihadrons/analysis.C b/macro/dihadrons/analysis.C new file mode 100644 index 00000000..d1d16139 --- /dev/null +++ b/macro/dihadrons/analysis.C @@ -0,0 +1,44 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Christopher Dilks + +R__LOAD_LIBRARY(EpicAnalysis) + +void analysis( + TString source="delphes", + Double_t eleBeamEn=18, Double_t ionBeamEn=275 + ) +{ + Analysis *A; + TString configFile, outfilePrefix; + if(source=="delphes") { + configFile = Form("datarec/hepmc.pythia6/radcor/%dx%d/files.config",(int)eleBeamEn,(int)ionBeamEn); + outfilePrefix = Form("dihadrons.delphes.%dx%d",(int)eleBeamEn,(int)ionBeamEn); + A = new AnalysisDelphes(configFile, outfilePrefix); + } else { + fmt::print(stderr,"ERROR: source '{}' not implemented\n",source); + return; + } + + // A->maxEvents = 30000; // limiter + A->writeSimpleTree = true; + A->SetReconMethod("Ele"); + A->includeOutputSet["1h"] = true; // optionally output single-hadron plots + + // dihadron final state ================================== + A->includeOutputSet["2h"] = true; // include the output set + A->AddFinalState("pipTrack"); // call `AddFinalState` exactly 2 times: once for each hadron + A->AddFinalState("pimTrack"); + + // cuts ================================================== + // - inclusive cuts + A->AddBinScheme("w"); A->BinScheme("w")->BuildBin("Min",3.0); // W > 3 GeV + A->AddBinScheme("y"); A->BinScheme("y")->BuildBin("Range",0.01,0.95); // 0.01 < y < 0.95 + // - single-hadron cuts: + /* don't use these for the dihadron final state, since they will only apply to the second hadron + A->AddBinScheme("z"); A->BinScheme("z")->BuildBin("Range",0.2,0.9); // 0.2 < z < 0.9 + A->AddBinScheme("xF"); A->BinScheme("xF")->BuildBin("Min",0.0); // xF > 0 + A->AddBinScheme("ptLab"); A->BinScheme("ptLab")->BuildBin("Min",0.1); // pT_lab > 0.1 GeV (tracking limit) + */ + + A->Execute(); +}; diff --git a/src/Analysis.cxx b/src/Analysis.cxx index 502e9804..091f42e4 100644 --- a/src/Analysis.cxx +++ b/src/Analysis.cxx @@ -53,6 +53,15 @@ Analysis::Analysis( availableBinSchemes.insert({ "JetEta", "jet eta" }); availableBinSchemes.insert({ "JetE", "jet energy" }); #endif + /* dihadrons */ + availableBinSchemes.insert({ "DihMh", "M_{h}" }); + availableBinSchemes.insert({ "DihMX", "M_{X}" }); + availableBinSchemes.insert({ "DihZ", "Z" }); + availableBinSchemes.insert({ "DihPhPerp", "P_{h,T}" }); + availableBinSchemes.insert({ "DihTheta", "#theta" }); + availableBinSchemes.insert({ "DihPhiH", "#phi_{h}" }); + availableBinSchemes.insert({ "DihPhiR", "#phi_{R}" }); + availableBinSchemes.insert({ "DihPhiS", "#phi_{S}" }); // available final states // - specify which final states you want to include using `AddFinalState(TString name)` @@ -87,6 +96,7 @@ Analysis::Analysis( // - the default settings are set here; override them at the macro level includeOutputSet.insert({ "inclusive", true }); // inclusive kinematics includeOutputSet.insert({ "1h", true }); // single hadron kinematics + includeOutputSet.insert({ "2h", false }); // dihadron kinematics includeOutputSet.insert({ "jets", false }); // jet kinematics includeOutputSet.insert({ "depolarization", false }); // depolarization factors & ratios @@ -107,6 +117,7 @@ Analysis::Analysis( weightInclusive = new WeightsUniform(); weightTrack = new WeightsUniform(); + weightDihadron = new WeightsUniform(); weightJet = new WeightsUniform(); // miscellaneous @@ -310,11 +321,11 @@ void Analysis::Prepare() { outFile = new TFile(outfileName,"RECREATE"); // instantiate shared objects - kin = new Kinematics(eleBeamEn,ionBeamEn,crossingAngle); + kin = new Kinematics(eleBeamEn, ionBeamEn, crossingAngle); kinTrue = new Kinematics(eleBeamEn, ionBeamEn, crossingAngle); - ST = new SimpleTree("tree",kin,kinTrue); - HFST = new HFSTree("hfstree",kin,kinTrue); - PT = new ParticleTree("ptree"); + ST = new SimpleTree("tree",kin,kinTrue); + HFST = new HFSTree("hfstree",kin,kinTrue); + PT = new ParticleTree("ptree"); // if including jets, define a `jet` final state #ifndef EXCLUDE_DELPHES @@ -328,9 +339,36 @@ void Analysis::Prepare() { includeOutputSet.insert({ "inclusive_only", includeOutputSet["inclusive"] && !includeOutputSet["1h"] + && !includeOutputSet["2h"] && !includeOutputSet["jets"] }); + // if including dihadrons, define a dihadron final state + if(includeOutputSet["2h"]) { + if(activeFinalStates.size()!=2) { + fmt::print(stderr,"ERROR: cannot include dihadron outputSet, since there should only be 2 final states defined\n"); + includeOutputSet["2h"] = false; + } else { + // add to dihSet + dihSet = new DihadronSet(); + // set finalStateID and title + TString dihadronFinalState = ""; + TString dihadronTitle = ""; + for(auto state : activeFinalStates) { + dihadronFinalState += state + "_"; + dihadronTitle += finalStateToTitle.at(state); + dihadronTitle(TRegexp(" .*")) = ""; + dihSet->IncludeHadron(state); + } + dihadronFinalState(TRegexp("_$")) = ""; + dihadronTitle += " dihadrons"; + // add the new dihadron final state + finalStateToTitle.insert({dihadronFinalState,dihadronTitle}); + AddFinalState(dihadronFinalState); + dihSet->SetFinalStateID(dihadronFinalState); + } + } + // if there are no final states defined, default to definitions here: if(BinScheme("finalState")->GetNumBins()==0) { std::cout << "NOTE: adding pi+ tracks for final state, since you specified none" << std::endl; @@ -369,6 +407,15 @@ void Analysis::Prepare() { HD->SetBinSchemeValue("phiS", [this](){ return kin->phiS; }); HD->SetBinSchemeValue("tSpin", [this](){ return (Double_t)kin->tSpin; }); HD->SetBinSchemeValue("lSpin", [this](){ return (Double_t)kin->lSpin; }); + /* dihadron */ + HD->SetBinSchemeValue("DihMh", [this](){ return dih->Mh; }); + HD->SetBinSchemeValue("DihMX", [this](){ return dih->MX; }); + HD->SetBinSchemeValue("DihZ", [this](){ return dih->Z; }); + HD->SetBinSchemeValue("DihPhPerp", [this](){ return dih->PhPerp; }); + HD->SetBinSchemeValue("DihTheta", [this](){ return dih->Theta; }); + HD->SetBinSchemeValue("DihPhiH", [this](){ return dih->PhiH; }); + HD->SetBinSchemeValue("DihPhiR", [this](){ return dih->PhiR; }); + HD->SetBinSchemeValue("DihPhiS", [this](){ return dih->PhiS; }); /* jets */ #ifndef EXCLUDE_DELPHES HD->SetBinSchemeValue("JetPT", [this](){ return kin->pTjet; }); @@ -499,7 +546,26 @@ void Analysis::Prepare() { NBINS,-TMath::Pi(),TMath::Pi() ); } - + // -- dihadron kinematics + if(includeOutputSet["2h"]) { + HS->DefineHist1D("DihMh", "M_{h}", "GeV", 2*NBINS, 0, 5); + HS->DefineHist1D("DihMX", "M_{X}", "GeV", NBINS, 0, 40); + HS->DefineHist1D("DihZ", "Z", "", NBINS, 0, 1); + HS->DefineHist1D("DihPhPerp", "P_{h,T}", "GeV", NBINS, 1e-2, 3, true); + HS->DefineHist1D("DihTheta", "#theta", "", NBINS, 0, TMath::Pi()); + HS->DefineHist1D("DihPhiH", "#phi_{h}", "", NBINS, -TMath::Pi(), TMath::Pi()); + HS->DefineHist1D("DihPhiR", "#phi_{R}", "", NBINS, -TMath::Pi(), TMath::Pi()); + HS->DefineHist1D("DihPhiS", "#phi_{S}", "", NBINS, -TMath::Pi(), TMath::Pi()); + HS->DefineHist2D("DihPhiHvsPhiR", "#phi_{R}", "#phi_{h}", "", "", + NBINS, -TMath::Pi(), TMath::Pi(), + NBINS, -TMath::Pi(), TMath::Pi() + ); + HS->DefineHist2D("DihThetaVsPh", "P_{h}", "#theta", "GeV", "", + NBINS, 1e-2, 40, + NBINS, 0, TMath::Pi(), + true + ); + } // -- jet kinematics #ifndef EXCLUDE_DELPHES if(includeOutputSet["jets"]) { @@ -562,6 +628,7 @@ void Analysis::Prepare() { // initialize total weights wInclusiveTotal = 0.; wTrackTotal = 0.; + wDihadronTotal = 0.; wJetTotal = 0.; }; @@ -670,10 +737,12 @@ void Analysis::Finish() { HD->Payload([this](Histos *H){ H->Write(); }); HD->ExecuteAndClearOps(); std::vector vec_wInclusiveTotal { wInclusiveTotal }; std::vector vec_wTrackTotal { wTrackTotal }; + std::vector vec_wDihadronTotal { wDihadronTotal }; std::vector vec_wJetTotal { wJetTotal }; outFile->WriteObject(&Q2xsecs, "XsTotal"); outFile->WriteObject(&vec_wInclusiveTotal, "WeightInclusiveTotal"); outFile->WriteObject(&vec_wTrackTotal, "WeightTrackTotal"); + outFile->WriteObject(&vec_wDihadronTotal, "WeightDihadronTotal"); outFile->WriteObject(&vec_wJetTotal, "WeightJetTotal"); // write binning schemes @@ -730,11 +799,18 @@ void Analysis::AddFinalState(TString finalStateN) { return; }; BinScheme("finalState")->BuildExternalBin(finalStateN,finalStateT); - activeFinalStates.insert(finalStateN); + activeFinalStates.push_back(finalStateN); fmt::print("AddFinalState: name='{}'\n title='{}'\n",finalStateN,finalStateT); }; +// check if this final state bin has been added +//---------------------------------------------- +Bool_t Analysis::IsFinalState(TString finalState) { + return std::find( activeFinalStates.begin(), activeFinalStates.end(), finalState) != activeFinalStates.end(); +}; + + // FillHistos methods: check bins and fill associated histograms // - checks which bins the track/jet/etc. falls in // - fills the histograms in the associated Histos objects @@ -842,6 +918,23 @@ void Analysis::FillHistos1h(Double_t wgt) { FillHistos(fill_payload); }; +// fill 2h (dihadron) histograms +void Analysis::FillHistos2h(Double_t wgt) { + auto fill_payload = [this,wgt] (Histos *H) { + H->FillHist1D("DihMh", dih->Mh, wgt); + H->FillHist1D("DihMX", dih->MX, wgt); + H->FillHist1D("DihZ", dih->Z, wgt); + H->FillHist1D("DihPhPerp", dih->PhPerp, wgt); + H->FillHist1D("DihTheta", dih->Theta, wgt); + H->FillHist1D("DihPhiH", dih->PhiH, wgt); + H->FillHist1D("DihPhiR", dih->PhiR, wgt); + H->FillHist1D("DihPhiS", dih->PhiS, wgt); + H->FillHist2D("DihPhiHvsPhiR", dih->PhiR, dih->PhiH, wgt); + H->FillHist2D("DihThetaVsPh", dih->Ph, dih->Theta, wgt); + }; + FillHistos(fill_payload); +} + // fill jet histograms void Analysis::FillHistosJets(Double_t wgt) { #ifndef EXCLUDE_DELPHES diff --git a/src/Analysis.h b/src/Analysis.h index e98b1832..4ab6aec4 100644 --- a/src/Analysis.h +++ b/src/Analysis.h @@ -32,12 +32,16 @@ #include "Histos.h" #include "HistosDAG.h" #include "Kinematics.h" +#include "Dihadrons.h" #include "SimpleTree.h" #include "HFSTree.h" #include "ParticleTree.h" #include "Weights.h" #include "CommonConstants.h" +class Dihadron; +class DihadronSet; + class Analysis : public TNamed { public: @@ -58,6 +62,7 @@ class Analysis : public TNamed // add a new final state bin void AddFinalState(TString finalStateN); + Bool_t IsFinalState(TString finalState); // check if it's been added // common settings Bool_t verbose; // if true, print a lot more information @@ -122,6 +127,7 @@ class Analysis : public TNamed // `FillHistos(weight)` methods: fill histograms void FillHistosInclusive(Double_t wgt); // inclusive kinematics void FillHistos1h(Double_t wgt); // single-hadron kinematics + void FillHistos2h(Double_t wgt); // dihadron kinematics void FillHistosJets(Double_t wgt); // jet kinematics // shared objects @@ -129,11 +135,14 @@ class Analysis : public TNamed HFSTree *HFST; ParticleTree *PT; Kinematics *kin, *kinTrue; + Dihadron *dih, *dihTrue; + DihadronSet *dihSet; HistosDAG *HD; Weights const* weightInclusive; Weights const* weightTrack; + Weights const* weightDihadron; Weights const* weightJet; - Double_t wInclusiveTotal, wTrackTotal, wJetTotal; + Double_t wInclusiveTotal, wTrackTotal, wDihadronTotal, wJetTotal; Long64_t entriesTot; Long64_t errorCnt; const TString sep = "--------------------------------------------"; @@ -172,7 +181,7 @@ class Analysis : public TNamed std::map reconMethodToTitle; std::map finalStateToTitle; std::map PIDtoFinalState; - std::set activeFinalStates; + std::vector activeFinalStates; // check if Q2 `val` is between `min` and `max`; if `max==0`, only `val>=min` is checked template bool InQ2Range(T val, T min, T max, bool ignoreZero=false) { @@ -202,6 +211,9 @@ class Analysis : public TNamed fmt::print("}}\n"); } + // expose protected members to `DihadronSet`, such as `kin` and `kinTrue` + friend class DihadronSet; + private: // fill histograms, according to `fill_payload` diff --git a/src/AnalysisAthena.cxx b/src/AnalysisAthena.cxx index c7c1c909..f92de2ec 100644 --- a/src/AnalysisAthena.cxx +++ b/src/AnalysisAthena.cxx @@ -233,7 +233,7 @@ void AnalysisAthena::Execute() // histograms; if not, proceed to next track auto kv = PIDtoFinalState.find(pid_); if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else continue; - if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) continue; + if(!IsFinalState(finalStateID)) continue; // calculate reconstructed hadron kinematics kin->vecHadron = part.vecPart; diff --git a/src/AnalysisDelphes.cxx b/src/AnalysisDelphes.cxx index cfff9534..afae7017 100644 --- a/src/AnalysisDelphes.cxx +++ b/src/AnalysisDelphes.cxx @@ -183,7 +183,7 @@ void AnalysisDelphes::Execute() { ); auto kv = PIDtoFinalState.find(pid); if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else continue; - if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) continue; + if(!IsFinalState(finalStateID)) continue; // get parent particle, to check if pion is from vector meson GenParticle *trkParticle = (GenParticle*)trk->Particle.GetObject(); @@ -228,12 +228,19 @@ void AnalysisDelphes::Execute() { // - `IsActiveEvent()` is only true if at least one bin gets filled for this track if( writeSimpleTree && HD->IsActiveEvent() ) ST->FillTree(wTrack); } + if(includeOutputSet["2h"]) dihSet->AddHadron(this); // tests //kin->ValidateHeadOnFrame(); }; // end track loop + // dihadrons - - - - - - - - - - - - - - - - - - - - - - - - - - - - + if(includeOutputSet["2h"]) { + auto wDihadron = Q2weightFactor * weightDihadron->GetWeight(*kinTrue); + wDihadronTotal += wDihadron; + dihSet->CalculateKinematics(this, wDihadron); // calls FillHistos2h(wDihadron) and FillHistosInclusive(wDihadron); + } // jet loop - - - - - - - - - - - - - - - - - - - - - - - - - - - - if(includeOutputSet["jets"]) { diff --git a/src/AnalysisEcce.cxx b/src/AnalysisEcce.cxx index ae78b325..efff1068 100644 --- a/src/AnalysisEcce.cxx +++ b/src/AnalysisEcce.cxx @@ -448,7 +448,7 @@ void AnalysisEcce::Execute() // histograms; if not, proceed to next track auto kv = PIDtoFinalState.find(pid_); if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else continue; - if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) continue; + if(!IsFinalState(finalStateID)) continue; // calculate reconstructed hadron kinematics kin->vecHadron = part.vecPart; diff --git a/src/AnalysisEpic.cxx b/src/AnalysisEpic.cxx index 47d8af5c..8341214d 100644 --- a/src/AnalysisEpic.cxx +++ b/src/AnalysisEpic.cxx @@ -322,7 +322,7 @@ void AnalysisEpic::Execute() // histograms; if not, proceed to next track auto kv = PIDtoFinalState.find(pid_); if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else continue; - if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) continue; + if(!IsFinalState(finalStateID)) continue; // calculate reconstructed hadron kinematics kin->vecHadron = part.vecPart; diff --git a/src/AnalysisEpicPodio.cxx b/src/AnalysisEpicPodio.cxx index 245fb185..1bcbdd0a 100644 --- a/src/AnalysisEpicPodio.cxx +++ b/src/AnalysisEpicPodio.cxx @@ -241,7 +241,7 @@ void AnalysisEpicPodio::Execute() // - check PID, to see if it's a final state we're interested in auto kv = PIDtoFinalState.find(recPDG); if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else return; - if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) return; + if(!IsFinalState(finalStateID)) return; // set SIDIS particle 4-momenta, and calculate their kinematics kinTrue->vecHadron = GetP4(simPart); diff --git a/src/Dihadrons.cxx b/src/Dihadrons.cxx new file mode 100644 index 00000000..3417759b --- /dev/null +++ b/src/Dihadrons.cxx @@ -0,0 +1,123 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Christopher Dilks + +#include "Dihadrons.h" + +ClassImp(DihadronSet) +ClassImp(Dihadron) + +Dihadron::Dihadron(TLorentzVector vecHad1, TLorentzVector vecHad2, Kinematics *K) +{ + const TLorentzVector &vecE = K->vecEleBeam; + const TLorentzVector &vecP = K->vecIonBeam; + const TLorentzVector &vecL = K->vecElectron; + const TLorentzVector &vecQ = K->vecQ; + const TLorentzVector &vecW = K->vecW; + const TLorentzVector &vecPh = vecHad1 + vecHad2; + const TLorentzVector &vecR = 0.5 * (vecHad1-vecHad2); + const TLorentzVector vecH[2] = {vecHad1, vecHad2}; + + // boost to ion rest frame + TLorentzVector IvecH[2]; + TLorentzVector IvecL, IvecQ, IvecPh; + K->BoostToIonFrame(vecL,IvecL); + K->BoostToIonFrame(vecQ,IvecQ); + K->BoostToIonFrame(vecPh,IvecPh); + for(int h=0; h<2; h++) + K->BoostToIonFrame(vecH[h],IvecH[h]); + + // invariant mass, missing mass, Z, PhPerp, Zeta + Mh = vecPh.M(); + MX = TMath::Abs((vecW-vecPh).M()); + Z = vecP.Dot(vecPh) / vecP.Dot(vecQ); + for(int h=0; h<2; h++) + hadZ[h] = vecP.Dot(vecH[h]) / vecP.Dot(vecQ); + Ph = vecPh.Vect().Mag(); + PhPerp = Kinematics::Reject(IvecPh.Vect(),IvecQ.Vect()).Mag(); + Zeta = 2 * vecR.Dot(vecP) / vecPh.Dot(vecP); + + // PhiH + PhiH = Kinematics::PlaneAngle( + IvecQ.Vect(), IvecL.Vect(), + IvecQ.Vect(), IvecPh.Vect() + ); + + // PhiR + Double_t coeff = + K->x * (Zeta*Mh*Mh - (vecH[0].M2() - vecH[1].M2()) ) / + ( K->Q2 * Z ); + TLorentzVector vecRT = vecR - (Zeta/2)*vecPh + coeff*vecP; + TLorentzVector IvecRT; + K->BoostToIonFrame(vecRT,IvecRT); + PhiR = Kinematics::PlaneAngle( + IvecQ.Vect(), IvecL.Vect(), + IvecQ.Vect(), IvecRT.Vect() + ); + + // PhiS + PhiS = K->phiS; + + // Theta + Double_t MRterm[2]; + for(int h=0; h<2; h++) + MRterm[h] = TMath::Sqrt( vecH[h].M2() + vecR.Vect().Mag2() ); + Theta = TMath::ACos( + ( MRterm[0] - MRterm[1] - Mh*Zeta ) / + ( 2 * vecR.Vect().Mag() ) + ); +} + +//////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////// + +void DihadronSet::IncludeHadron(TString hadName) { + includedHadrons.push_back(hadName); +} + +void DihadronSet::AddHadron(Analysis *A) { + if(hadSetRec.find(A->finalStateID)==hadSetRec.end()) { + hadSetRec.insert({A->finalStateID,{}}); + hadSetGen.insert({A->finalStateID,{}}); + } + if(debug) fmt::print("DihadronSet: AddHadron '{}'\n",A->finalStateID); + hadSetRec[A->finalStateID].push_back( A->kin->vecHadron ); + hadSetGen[A->finalStateID].push_back( A->kinTrue->vecHadron ); +} + +void DihadronSet::CalculateKinematics(Analysis *A, Double_t wgt) { + + // TODO: this does not yet handle generalized dihadrons + if(includedHadrons.size()!=2) fmt::print("ERROR: more or less than 2 final states defined for DihadronSet\n"); + + // hadron pairing + auto PairHadrons = [this] (auto hadSet, auto &dihList, Kinematics *K) { + TString hadNames[2] = { includedHadrons[0], includedHadrons[1] }; + for(const auto &had0vec : hadSet[hadNames[0]]) { + for(const auto &had1vec : hadSet[hadNames[1]]) { + dihList.push_back(Dihadron(had0vec, had1vec, K)); + if(debug) fmt::print(" pair: {}, {}\n",hadNames[0],hadNames[1]); + } + } + }; + if(debug) fmt::print("DihadronSet::CalculateKinematics reconstructed\n"); + PairHadrons(hadSetRec,dihListRec,A->kin); + if(debug) fmt::print("DihadronSet::CalculateKinematics generated\n"); + PairHadrons(hadSetGen,dihListGen,A->kinTrue); + + // fill output data structures + for(std::size_t i=0; idih = &(dihListRec[i]); + A->dihTrue = &(dihListGen[i]); + auto finalStateID_tmp = A->finalStateID; // temporarily change Analysis::finalStateID + A->finalStateID = dihadronFinalStateID; // to that of this DihadronSet + A->FillHistos2h(wgt); + A->FillHistosInclusive(wgt); + A->finalStateID = finalStateID_tmp; // revert Analysis::finalStateID + } + + // reset internal storage + hadSetRec.clear(); + hadSetGen.clear(); + dihListRec.clear(); + dihListGen.clear(); +} diff --git a/src/Dihadrons.h b/src/Dihadrons.h new file mode 100644 index 00000000..fb483975 --- /dev/null +++ b/src/Dihadrons.h @@ -0,0 +1,39 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Christopher Dilks + +#pragma once + +#include +#include "Analysis.h" +#include "Kinematics.h" + +class Analysis; + +class Dihadron { + public: + Dihadron() {}; + Dihadron(TLorentzVector vecHad1, TLorentzVector vecHad2, Kinematics *K); + ~Dihadron() {}; + Double_t Mh, MX, Z, Ph, PhPerp, Theta, PhiH, PhiR, PhiS, Zeta; + Double_t hadZ[2]; + private: + ClassDef(Dihadron,1); +}; + +class DihadronSet { + public: + DihadronSet() : debug(false), dihadronFinalStateID("") {}; + ~DihadronSet() {}; + void IncludeHadron(TString hadName); + void AddHadron(Analysis *A); + void SetFinalStateID(TString state) { dihadronFinalStateID = state; } + void CalculateKinematics(Analysis *A, Double_t wgt=1.0); + private: + std::vector includedHadrons; + std::map> hadSetRec, hadSetGen; + std::vector dihListRec, dihListGen; + TString dihadronFinalStateID; + bool debug; + ClassDef(DihadronSet,1); +}; + diff --git a/src/Kinematics.h b/src/Kinematics.h index aee5e953..a5b80d6f 100644 --- a/src/Kinematics.h +++ b/src/Kinematics.h @@ -204,6 +204,12 @@ class Kinematics : public TObject // - convert energy,mass to momentum static Double_t EMtoP(Double_t energy, Double_t mass) { return TMath::Sqrt( TMath::Power(energy,2) - TMath::Power(mass,2) ); + }; + // - get angle between two vectors + static Double_t AngleSubtend(TVector3 vA, TVector3 vB) { + Double_t m = vA.Mag() * vB.Mag(); + if(m>0) return TMath::ACos( vA.Dot(vB) / m ); + return -10000; }; // - vector projection: returns vA projected onto vB static TVector3 Project(TVector3 vA, TVector3 vB) { diff --git a/src/LinkDef.h b/src/LinkDef.h index e9c023ff..49d23908 100644 --- a/src/LinkDef.h +++ b/src/LinkDef.h @@ -16,6 +16,8 @@ // analysis objects #pragma link C++ class Kinematics+; +#pragma link C++ class Dihadron+; +#pragma link C++ class DihadronSet+; #pragma link C++ class SimpleTree+; #pragma link C++ class HFSTree+; #pragma link C++ class ParticleTree+; diff --git a/tutorial/analysis_template.C b/tutorial/analysis_template.C index 35f82185..bcb9e740 100644 --- a/tutorial/analysis_template.C +++ b/tutorial/analysis_template.C @@ -32,6 +32,7 @@ void analysis_template( A->includeOutputSet["jets"] = true; // - additional example settings; see `src/Analysis.cxx` for more //A->includeOutputSet["1h"] = false; + //A->includeOutputSet["2h"] = true; //A->includeOutputSet["inclusive"] = false; //A->includeOutputSet["depolarization"] = true;