Skip to content

Commit 5d41c1c

Browse files
c-dilksc-dilks
committed
add most dihadron kinematics calculations
1 parent 4de779a commit 5d41c1c

File tree

10 files changed

+125
-27
lines changed

10 files changed

+125
-27
lines changed

‎macro/dihadrons/analysis.C‎

Lines changed: 9 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -20,18 +20,23 @@ void analysis(
2020
A->maxEvents = 30000; // limiter
2121
A->writeSimpleTree = true;
2222
A->SetReconMethod("Ele");
23-
A->includeOutputSet["1h"] = true;
24-
A->includeOutputSet["2h"] = true;
23+
A->includeOutputSet["1h"] = true; // optionally output single-hadron plots
2524

26-
A->AddFinalState("pipTrack");
25+
// dihadron final state ==================================
26+
A->includeOutputSet["2h"] = true; // include the output set
27+
A->AddFinalState("pipTrack"); // call `AddFinalState` exactly 2 times: once for each hadron
2728
A->AddFinalState("pimTrack");
2829

29-
// single-hadron cuts ====================================
30+
// cuts ==================================================
31+
// - inclusive cuts
3032
A->AddBinScheme("w"); A->BinScheme("w")->BuildBin("Min",3.0); // W > 3 GeV
3133
A->AddBinScheme("y"); A->BinScheme("y")->BuildBin("Range",0.01,0.95); // 0.01 < y < 0.95
34+
// - single-hadron cuts:
35+
/* don't use these for the dihadron final state, since they will only apply to the second hadron
3236
A->AddBinScheme("z"); A->BinScheme("z")->BuildBin("Range",0.2,0.9); // 0.2 < z < 0.9
3337
A->AddBinScheme("xF"); A->BinScheme("xF")->BuildBin("Min",0.0); // xF > 0
3438
A->AddBinScheme("ptLab"); A->BinScheme("ptLab")->BuildBin("Min",0.1); // pT_lab > 0.1 GeV (tracking limit)
39+
*/
3540

3641
A->Execute();
3742
};

‎src/Analysis.cxx‎

Lines changed: 41 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -68,7 +68,14 @@ Analysis::Analysis(
6868
PIDtoFinalState.insert({ 2212, "pTrack" });
6969

7070
// dihadrons
71-
availableBinSchemes.insert({ "dihMh", "M_{h}" });
71+
availableBinSchemes.insert({ "DihMh", "M_{h}" });
72+
availableBinSchemes.insert({ "DihMX", "M_{X}" });
73+
availableBinSchemes.insert({ "DihZ", "Z" });
74+
availableBinSchemes.insert({ "DihPhPerp", "P_{h,T}" });
75+
availableBinSchemes.insert({ "DihTheta", "#theta" });
76+
availableBinSchemes.insert({ "DihPhiH", "#phi_{h}" });
77+
availableBinSchemes.insert({ "DihPhiR", "#phi_{R}" });
78+
availableBinSchemes.insert({ "DihPhiS", "#phi_{S}" });
7279

7380
// kinematics reconstruction methods
7481
// - choose one of these methods using `SetReconMethod(TString name)`
@@ -340,11 +347,6 @@ void Analysis::Prepare() {
340347
dihSet->IncludeHadron(state);
341348
}
342349
dihadronFinalState(TRegexp("_$")) = "";
343-
// aesthetic quick fix: re-order the name and title for pi+pi-
344-
if(dihadronFinalState=="pimTrack_pipTrack") {
345-
dihadronFinalState = "pipTrack_pimTrack";
346-
dihadronTitle = "#pi^{+}#pi^{-}";
347-
}
348350
dihadronTitle += " dihadrons";
349351
// add the new dihadron final state
350352
finalStateToTitle.insert({dihadronFinalState,dihadronTitle});
@@ -381,7 +383,14 @@ void Analysis::Prepare() {
381383
HD->SetBinSchemeValue("tSpin", [this](){ return (Double_t)kin->tSpin; });
382384
HD->SetBinSchemeValue("lSpin", [this](){ return (Double_t)kin->lSpin; });
383385
/* dihadron */
384-
HD->SetBinSchemeValue("dihMh", [this](){ return dih->Mh; });
386+
HD->SetBinSchemeValue("DihMh", [this](){ return dih->Mh; });
387+
HD->SetBinSchemeValue("DihMX", [this](){ return dih->MX; });
388+
HD->SetBinSchemeValue("DihZ", [this](){ return dih->Z; });
389+
HD->SetBinSchemeValue("DihPhPerp", [this](){ return dih->PhPerp; });
390+
HD->SetBinSchemeValue("DihTheta", [this](){ return dih->Theta; });
391+
HD->SetBinSchemeValue("DihPhiH", [this](){ return dih->PhiH; });
392+
HD->SetBinSchemeValue("DihPhiR", [this](){ return dih->PhiR; });
393+
HD->SetBinSchemeValue("DihPhiS", [this](){ return dih->PhiS; });
385394
/* jets */
386395
#ifndef EXCLUDE_DELPHES
387396
HD->SetBinSchemeValue("JetPT", [this](){ return kin->pTjet; });
@@ -512,7 +521,14 @@ void Analysis::Prepare() {
512521
}
513522
// -- dihadron kinematics
514523
if(includeOutputSet["2h"]) {
515-
HS->DefineHist1D("dihMh", "M_{h}", "GeV", 2*NBINS, 0, 5);
524+
HS->DefineHist1D("DihMh", "M_{h}", "GeV", 2*NBINS, 0, 5);
525+
HS->DefineHist1D("DihMX", "M_{X}", "GeV", NBINS, 0, 40);
526+
HS->DefineHist1D("DihZ", "Z", "", NBINS, 0, 1);
527+
HS->DefineHist1D("DihPhPerp", "P_{h,T}", "GeV", NBINS, 1e-2, 3, true);
528+
HS->DefineHist1D("DihTheta", "#theta", "", NBINS, 0, TMath::Pi());
529+
HS->DefineHist1D("DihPhiH", "#phi_{h}", "", NBINS, -TMath::Pi(), TMath::Pi());
530+
HS->DefineHist1D("DihPhiR", "#phi_{R}", "", NBINS, -TMath::Pi(), TMath::Pi());
531+
HS->DefineHist1D("DihPhiS", "#phi_{S}", "", NBINS, -TMath::Pi(), TMath::Pi());
516532
}
517533
// -- jet kinematics
518534
#ifndef EXCLUDE_DELPHES
@@ -723,11 +739,19 @@ void Analysis::AddFinalState(TString finalStateN) {
723739
return;
724740
};
725741
BinScheme("finalState")->BuildExternalBin(finalStateN,finalStateT);
726-
activeFinalStates.insert(finalStateN);
742+
activeFinalStates.push_back(finalStateN);
727743
fmt::print("AddFinalState: name='{}'\n title='{}'\n",finalStateN,finalStateT);
728744
};
729745

730746

747+
// check if this final state bin has been added
748+
//----------------------------------------------
749+
Bool_t Analysis::IsFinalState(TString finalState) {
750+
return std::find( activeFinalStates.begin(), activeFinalStates.end(), finalState) != activeFinalStates.end();
751+
>>>>>>> add most dihadron kinematics calculations
752+
};
753+
754+
731755
// FillHistos methods: check bins and fill associated histograms
732756
// - checks which bins the track/jet/etc. falls in
733757
// - fills the histograms in the associated Histos objects
@@ -828,7 +852,14 @@ void Analysis::FillHistosTracks() {
828852
void Analysis::FillHistosDihadrons() {
829853
HD->CheckBins();
830854
HD->Payload([this](Histos *H){
831-
H->FillHist1D("dihMh", dih->Mh, wTrack);
855+
H->FillHist1D("DihMh", dih->Mh, wTrack);
856+
H->FillHist1D("DihMX", dih->MX, wTrack);
857+
H->FillHist1D("DihZ", dih->Z, wTrack);
858+
H->FillHist1D("DihPhPerp", dih->PhPerp, wTrack);
859+
H->FillHist1D("DihTheta", dih->Theta, wTrack);
860+
H->FillHist1D("DihPhiH", dih->PhiH, wTrack);
861+
H->FillHist1D("DihPhiR", dih->PhiR, wTrack);
862+
H->FillHist1D("DihPhiS", dih->PhiS, wTrack);
832863
});
833864
HD->ExecuteOps(true);
834865
}

‎src/Analysis.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -56,6 +56,7 @@ class Analysis : public TNamed
5656

5757
// add a new final state bin
5858
void AddFinalState(TString finalStateN);
59+
Bool_t IsFinalState(TString finalState); // check if it's been added
5960

6061
// common settings
6162
Bool_t verbose; // if true, print a lot more information
@@ -165,7 +166,7 @@ class Analysis : public TNamed
165166
std::map<TString,TString> reconMethodToTitle;
166167
std::map<TString,TString> finalStateToTitle;
167168
std::map<int,TString> PIDtoFinalState;
168-
std::set<TString> activeFinalStates;
169+
std::vector<TString> activeFinalStates;
169170

170171
// check if Q2 `val` is between `min` and `max`; if `max==0`, only `val>=min` is checked
171172
template<class T> bool InQ2Range(T val, T min, T max, bool ignoreZero=false) {

‎src/AnalysisAthena.cxx‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -220,7 +220,7 @@ void AnalysisAthena::Execute()
220220
// histograms; if not, proceed to next track
221221
auto kv = PIDtoFinalState.find(pid_);
222222
if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else continue;
223-
if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) continue;
223+
if(!IsFinalState(finalStateID)) continue;
224224

225225
// calculate reconstructed hadron kinematics
226226
kin->vecHadron = part.vecPart;

‎src/AnalysisDelphes.cxx‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -174,7 +174,7 @@ void AnalysisDelphes::Execute() {
174174
);
175175
auto kv = PIDtoFinalState.find(pid);
176176
if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else continue;
177-
if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) continue;
177+
if(!IsFinalState(finalStateID)) continue;
178178

179179
// get parent particle, to check if pion is from vector meson
180180
GenParticle *trkParticle = (GenParticle*)trk->Particle.GetObject();

‎src/AnalysisEcce.cxx‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -306,7 +306,7 @@ void AnalysisEcce::Execute()
306306
// histograms; if not, proceed to next track
307307
auto kv = PIDtoFinalState.find(pid_);
308308
if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else continue;
309-
if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) continue;
309+
if(!IsFinalState(finalStateID)) continue;
310310

311311
// calculate reconstructed hadron kinematics
312312
kin->vecHadron = part.vecPart;

‎src/AnalysisEpic.cxx‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -225,7 +225,7 @@ void AnalysisEpic::Execute()
225225
// - check PID, to see if it's a final state we're interested in
226226
auto kv = PIDtoFinalState.find(recPDG);
227227
if(kv!=PIDtoFinalState.end()) finalStateID = kv->second; else return;
228-
if(activeFinalStates.find(finalStateID)==activeFinalStates.end()) return;
228+
if(!IsFinalState(finalStateID)) return;
229229

230230
// set SIDIS particle 4-momenta, and calculate their kinematics
231231
kinTrue->vecHadron = GetP4(simPart);

‎src/Dihadrons.cxx‎

Lines changed: 59 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -3,9 +3,62 @@
33
ClassImp(DihadronSet)
44
ClassImp(Dihadron)
55

6-
Dihadron::Dihadron(TLorentzVector vecHad1, TLorentzVector vecHad2) {
7-
auto vecPh = vecHad1 + vecHad2;
6+
Dihadron::Dihadron(TLorentzVector vecHad1, TLorentzVector vecHad2, Kinematics *K)
7+
{
8+
const TLorentzVector &vecEB = K->vecEleBeam;
9+
const TLorentzVector &vecIB = K->vecIonBeam;
10+
const TLorentzVector &vecL = K->vecElectron;
11+
const TLorentzVector &vecQ = K->vecQ;
12+
const TLorentzVector &vecW = K->vecW;
13+
const TLorentzVector &vecPh = vecHad1 + vecHad2;
14+
const TLorentzVector &vecR = 0.5 * (vecHad1-vecHad2);
15+
const TLorentzVector vecH[2] = {vecHad1, vecHad2};
16+
17+
// boost to ion rest frame
18+
TLorentzVector IvecH[2];
19+
TLorentzVector IvecL, IvecQ, IvecPh;
20+
K->BoostToIonFrame(vecL,IvecL);
21+
K->BoostToIonFrame(vecQ,IvecQ);
22+
K->BoostToIonFrame(vecPh,IvecPh);
23+
for(int h=0; h<2; h++)
24+
K->BoostToIonFrame(vecH[h],IvecH[h]);
25+
26+
// invariant mass, missing mass, Z, PhPerp
827
Mh = vecPh.M();
28+
MX = TMath::Abs((vecW-vecPh).M());
29+
Z = vecIB.Dot(vecPh) / vecIB.Dot(vecQ);
30+
for(int h=0; h<2; h++)
31+
hadZ[h] = vecIB.Dot(vecH[h]) / vecIB.Dot(vecQ);
32+
PhPerp = Kinematics::Reject(IvecPh.Vect(),IvecQ.Vect()).Mag();
33+
34+
// Theta
35+
auto boost_com = -1 * vecPh.BoostVector(); // boost to dihadron CoM frame
36+
TLorentzVector vecH_com[2];
37+
for(int h=0; h<2; h++) {
38+
vecH_com[h] = vecH[h];
39+
vecH_com[h].Boost(boost_com);
40+
}
41+
Theta = Kinematics::AngleSubtend(vecH_com[0].Vect(),vecPh.Vect());
42+
43+
// PhiH
44+
PhiH = Kinematics::PlaneAngle(
45+
vecQ.Vect(), vecL.Vect(),
46+
vecQ.Vect(), vecPh.Vect()
47+
);
48+
49+
// PhiR (calculated in ion rest frame)
50+
TVector3 momHad_perp[2]; // perp-frame hadron 3-momenta
51+
for(int h=0; h<2; h++)
52+
momHad_perp[h] = Kinematics::Reject(IvecH[h].Vect(), IvecQ.Vect());
53+
auto momRT = 1 / ( hadZ[0] + hadZ[1] ) *
54+
( hadZ[1] * momHad_perp[0] - hadZ[0] * momHad_perp[1] );
55+
PhiR = Kinematics::PlaneAngle(
56+
IvecQ.Vect(), IvecL.Vect(),
57+
IvecQ.Vect(), momRT
58+
);
59+
60+
// PhiS
61+
PhiS = K->phiS;
962
}
1063

1164
////////////////////////////////////////////////////////////////////
@@ -31,19 +84,19 @@ void DihadronSet::CalculateKinematics(Analysis *A) {
3184
if(includedHadrons.size()!=2) fmt::print("ERROR: more or less than 2 final states defined for DihadronSet\n");
3285

3386
// hadron pairing
34-
auto PairHadrons = [this] (auto hadSet, auto &dihList) {
87+
auto PairHadrons = [this] (auto hadSet, auto &dihList, Kinematics *K) {
3588
TString hadNames[2] = { includedHadrons[0], includedHadrons[1] };
3689
for(const auto &had0vec : hadSet[hadNames[0]]) {
3790
for(const auto &had1vec : hadSet[hadNames[1]]) {
38-
dihList.push_back(Dihadron(had0vec,had1vec));
91+
dihList.push_back(Dihadron(had0vec, had1vec, K));
3992
if(debug) fmt::print(" pair: {}, {}\n",hadNames[0],hadNames[1]);
4093
}
4194
}
4295
};
4396
if(debug) fmt::print("DihadronSet::CalculateKinematics reconstructed\n");
44-
PairHadrons(hadSetRec,dihListRec);
97+
PairHadrons(hadSetRec,dihListRec,A->kin);
4598
if(debug) fmt::print("DihadronSet::CalculateKinematics generated\n");
46-
PairHadrons(hadSetGen,dihListGen);
99+
PairHadrons(hadSetGen,dihListGen,A->kinTrue);
47100

48101
// fill output data structures
49102
for(std::size_t i=0; i<dihListRec.size(); i++) {

‎src/Dihadrons.h‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2,15 +2,17 @@
22

33
#include <vector>
44
#include "Analysis.h"
5+
#include "Kinematics.h"
56

67
class Analysis;
78

89
class Dihadron {
910
public:
1011
Dihadron() {};
11-
Dihadron(TLorentzVector vecHad1, TLorentzVector vecHad2);
12+
Dihadron(TLorentzVector vecHad1, TLorentzVector vecHad2, Kinematics *K);
1213
~Dihadron() {};
13-
Double_t Mh;
14+
Double_t Mh, MX, Z, PhPerp, Theta, PhiH, PhiR, PhiS;
15+
Double_t hadZ[2];
1416
private:
1517
ClassDef(Dihadron,1);
1618
};

‎src/Kinematics.h‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -169,6 +169,12 @@ class Kinematics : public TObject
169169
// - convert energy,mass to momentum
170170
static Double_t EMtoP(Double_t energy, Double_t mass) {
171171
return TMath::Sqrt( TMath::Power(energy,2) - TMath::Power(mass,2) );
172+
};
173+
// - get angle between two vectors
174+
static Double_t AngleSubtend(TVector3 vA, TVector3 vB) {
175+
Double_t m = vA.Mag() * vB.Mag();
176+
if(m>0) return TMath::ACos( vA.Dot(vB) / m );
177+
return -10000;
172178
};
173179
// - vector projection: returns vA projected onto vB
174180
static TVector3 Project(TVector3 vA, TVector3 vB) {

0 commit comments

Comments
 (0)