Skip to content

Commit

Permalink
add most dihadron kinematics calculations
Browse files Browse the repository at this point in the history
  • Loading branch information
@c-dilks
c-dilks committed Oct 22, 2022
1 parent aabbda0 commit bb01842
Show file tree
Hide file tree
Showing 10 changed files with 128 additions and 33 deletions.
13 changes: 9 additions & 4 deletions macro/dihadrons/analysis.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();
};
54 changes: 42 additions & 12 deletions src/Analysis.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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});
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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();
};


Expand Down Expand Up @@ -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);
}
Expand Down
3 changes: 2 additions & 1 deletion src/Analysis.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -165,7 +166,7 @@ class Analysis : public TNamed
std::map<TString,TString> reconMethodToTitle;
std::map<TString,TString> finalStateToTitle;
std::map<int,TString> PIDtoFinalState;
std::set<TString> activeFinalStates;
std::vector<TString> activeFinalStates;

// check if Q2 `val` is between `min` and `max`; if `max==0`, only `val>=min` is checked
template<class T> bool InQ2Range(T val, T min, T max, bool ignoreZero=false) {
Expand Down
2 changes: 1 addition & 1 deletion src/AnalysisAthena.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
8 changes: 3 additions & 5 deletions src/AnalysisDelphes.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down Expand Up @@ -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);
Expand Down
2 changes: 1 addition & 1 deletion src/AnalysisEcce.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
2 changes: 1 addition & 1 deletion src/AnalysisEpic.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
65 changes: 59 additions & 6 deletions src/Dihadrons.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}

////////////////////////////////////////////////////////////////////
Expand All @@ -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<dihListRec.size(); i++) {
Expand Down
6 changes: 4 additions & 2 deletions src/Dihadrons.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,15 +2,17 @@

#include <vector>
#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);
};
Expand Down
6 changes: 6 additions & 0 deletions src/Kinematics.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down

0 comments on commit bb01842

Please sign in to comment.