From bb0184266362f4c74ebc5a40fe5268ea6a3a9187 Mon Sep 17 00:00:00 2001 From: Christopher Dilks Date: Sat, 22 Oct 2022 02:33:47 -0400 Subject: [PATCH] add most dihadron kinematics calculations --- macro/dihadrons/analysis.C | 13 +++++--- src/Analysis.cxx | 54 ++++++++++++++++++++++++------- src/Analysis.h | 3 +- src/AnalysisAthena.cxx | 2 +- src/AnalysisDelphes.cxx | 8 ++--- src/AnalysisEcce.cxx | 2 +- src/AnalysisEpic.cxx | 2 +- src/Dihadrons.cxx | 65 ++++++++++++++++++++++++++++++++++---- src/Dihadrons.h | 6 ++-- src/Kinematics.h | 6 ++++ 10 files changed, 128 insertions(+), 33 deletions(-) diff --git a/macro/dihadrons/analysis.C b/macro/dihadrons/analysis.C index 01d5dbd3..41b6b897 100644 --- a/macro/dihadrons/analysis.C +++ b/macro/dihadrons/analysis.C @@ -20,18 +20,23 @@ void analysis( A->maxEvents = 30000; // limiter A->writeSimpleTree = true; A->SetReconMethod("Ele"); - A->includeOutputSet["1h"] = true; - A->includeOutputSet["2h"] = true; + A->includeOutputSet["1h"] = true; // optionally output single-hadron plots - A->AddFinalState("pipTrack"); + // 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"); - // single-hadron cuts ==================================== + // 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 b09ebb56..e957ef72 100644 --- a/src/Analysis.cxx +++ b/src/Analysis.cxx @@ -64,7 +64,14 @@ Analysis::Analysis( PIDtoFinalState.insert({ 2212, "pTrack" }); // dihadrons - availableBinSchemes.insert({ "dihMh", "M_{h}" }); + 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}" }); // jets #ifndef EXCLUDE_DELPHES @@ -338,11 +345,6 @@ void Analysis::Prepare() { dihSet->IncludeHadron(state); } dihadronFinalState(TRegexp("_$")) = ""; - // aesthetic quick fix: re-order the name and title for pi+pi- - if(dihadronFinalState=="pimTrack_pipTrack") { - dihadronFinalState = "pipTrack_pimTrack"; - dihadronTitle = "#pi^{+}#pi^{-}"; - } dihadronTitle += " dihadrons"; // add the new dihadron final state finalStateToTitle.insert({dihadronFinalState,dihadronTitle}); @@ -379,11 +381,18 @@ void Analysis::Prepare() { 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("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("ptJet", [this](){ return kin->pTjet; }); - HD->SetBinSchemeValue("zJet", [this](){ return kin->zjet; }); + HD->SetBinSchemeValue("ptJet", [this](){ return kin->pTjet; }); + HD->SetBinSchemeValue("zJet", [this](){ return kin->zjet; }); #endif // some bin schemes values are checked here, instead of by CutDef checks ("External" cut type) @@ -510,7 +519,14 @@ void Analysis::Prepare() { } // -- dihadron kinematics if(includeOutputSet["2h"]) { - HS->DefineHist1D("dihMh", "M_{h}", "GeV", 2*NBINS, 0, 5); + 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()); } // -- jet kinematics #ifndef EXCLUDE_DELPHES @@ -721,7 +737,14 @@ void Analysis::AddFinalState(TString finalStateN) { return; }; BinScheme("finalState")->BuildExternalBin(finalStateN,finalStateT); - activeFinalStates.insert(finalStateN); + activeFinalStates.push_back(finalStateN); +}; + + +// 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(); }; @@ -825,7 +848,14 @@ void Analysis::FillHistosTracks() { void Analysis::FillHistosDihadrons() { HD->CheckBins(); HD->Payload([this](Histos *H){ - H->FillHist1D("dihMh", dih->Mh, wTrack); + H->FillHist1D("DihMh", dih->Mh, wTrack); + H->FillHist1D("DihMX", dih->MX, wTrack); + H->FillHist1D("DihZ", dih->Z, wTrack); + H->FillHist1D("DihPhPerp", dih->PhPerp, wTrack); + H->FillHist1D("DihTheta", dih->Theta, wTrack); + H->FillHist1D("DihPhiH", dih->PhiH, wTrack); + H->FillHist1D("DihPhiR", dih->PhiR, wTrack); + H->FillHist1D("DihPhiS", dih->PhiS, wTrack); }); HD->ExecuteOps(true); } diff --git a/src/Analysis.h b/src/Analysis.h index 21e12f4e..7d957bd7 100644 --- a/src/Analysis.h +++ b/src/Analysis.h @@ -56,6 +56,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 @@ -165,7 +166,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) { diff --git a/src/AnalysisAthena.cxx b/src/AnalysisAthena.cxx index 969746e1..8ef1effb 100644 --- a/src/AnalysisAthena.cxx +++ b/src/AnalysisAthena.cxx @@ -220,7 +220,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 404bc22f..d14fa6fc 100644 --- a/src/AnalysisDelphes.cxx +++ b/src/AnalysisDelphes.cxx @@ -173,7 +173,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(); @@ -231,10 +231,8 @@ void AnalysisDelphes::Execute() { if(includeOutputSet["2h"]) dihSet->CalculateKinematics(this); // jet loop - - - - - - - - - - - - - - - - - - - - - - - - - - - - - finalStateID = "jet"; - // FIXME: probably don't need both a `finalStateID` and `includeOutputSet` to control whether - // we do anything with jets or not - if(activeFinalStates.find(finalStateID)!=activeFinalStates.end() && includeOutputSet["jets"]) { + if(includeOutputSet["jets"]) { + finalStateID = "jet"; #ifdef INCLUDE_CENTAURO if(useBreitJets) kin->GetBreitFrameJets(itEFlowTrack, itEFlowPhoton, itEFlowNeutralHadron, itParticle); diff --git a/src/AnalysisEcce.cxx b/src/AnalysisEcce.cxx index 956a3cf0..8475f162 100644 --- a/src/AnalysisEcce.cxx +++ b/src/AnalysisEcce.cxx @@ -306,7 +306,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 ddebb6d3..aaf75d9b 100644 --- a/src/AnalysisEpic.cxx +++ b/src/AnalysisEpic.cxx @@ -225,7 +225,7 @@ void AnalysisEpic::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 index de0d4b7c..cee06f83 100644 --- a/src/Dihadrons.cxx +++ b/src/Dihadrons.cxx @@ -3,9 +3,62 @@ ClassImp(DihadronSet) ClassImp(Dihadron) -Dihadron::Dihadron(TLorentzVector vecHad1, TLorentzVector vecHad2) { - auto vecPh = vecHad1 + vecHad2; +Dihadron::Dihadron(TLorentzVector vecHad1, TLorentzVector vecHad2, Kinematics *K) +{ + const TLorentzVector &vecEB = K->vecEleBeam; + const TLorentzVector &vecIB = 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 Mh = vecPh.M(); + MX = TMath::Abs((vecW-vecPh).M()); + Z = vecIB.Dot(vecPh) / vecIB.Dot(vecQ); + for(int h=0; h<2; h++) + hadZ[h] = vecIB.Dot(vecH[h]) / vecIB.Dot(vecQ); + PhPerp = Kinematics::Reject(IvecPh.Vect(),IvecQ.Vect()).Mag(); + + // Theta + auto boost_com = -1 * vecPh.BoostVector(); // boost to dihadron CoM frame + TLorentzVector vecH_com[2]; + for(int h=0; h<2; h++) { + vecH_com[h] = vecH[h]; + vecH_com[h].Boost(boost_com); + } + Theta = Kinematics::AngleSubtend(vecH_com[0].Vect(),vecPh.Vect()); + + // PhiH + PhiH = Kinematics::PlaneAngle( + vecQ.Vect(), vecL.Vect(), + vecQ.Vect(), vecPh.Vect() + ); + + // PhiR (calculated in ion rest frame) + TVector3 momHad_perp[2]; // perp-frame hadron 3-momenta + for(int h=0; h<2; h++) + momHad_perp[h] = Kinematics::Reject(IvecH[h].Vect(), IvecQ.Vect()); + auto momRT = 1 / ( hadZ[0] + hadZ[1] ) * + ( hadZ[1] * momHad_perp[0] - hadZ[0] * momHad_perp[1] ); + PhiR = Kinematics::PlaneAngle( + IvecQ.Vect(), IvecL.Vect(), + IvecQ.Vect(), momRT + ); + + // PhiS + PhiS = K->phiS; } //////////////////////////////////////////////////////////////////// @@ -31,19 +84,19 @@ void DihadronSet::CalculateKinematics(Analysis *A) { 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) { + 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)); + 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); + PairHadrons(hadSetRec,dihListRec,A->kin); if(debug) fmt::print("DihadronSet::CalculateKinematics generated\n"); - PairHadrons(hadSetGen,dihListGen); + PairHadrons(hadSetGen,dihListGen,A->kinTrue); // fill output data structures for(std::size_t i=0; i #include "Analysis.h" +#include "Kinematics.h" class Analysis; class Dihadron { public: Dihadron() {}; - Dihadron(TLorentzVector vecHad1, TLorentzVector vecHad2); + Dihadron(TLorentzVector vecHad1, TLorentzVector vecHad2, Kinematics *K); ~Dihadron() {}; - Double_t Mh; + Double_t Mh, MX, Z, PhPerp, Theta, PhiH, PhiR, PhiS; + Double_t hadZ[2]; private: ClassDef(Dihadron,1); }; diff --git a/src/Kinematics.h b/src/Kinematics.h index 431ac2bb..1c0d6a54 100644 --- a/src/Kinematics.h +++ b/src/Kinematics.h @@ -169,6 +169,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) {