libraries/ANALYSIS/DHistogramActions

Bug Summary

File:	/volatile/halld/gluex/nightly/2024-07-31/Linux_CentOS7.7-x86_64-gcc4.8.5/halld_recon/src/libraries/ANALYSIS/DHistogramActions_Reaction.cc
Location:	line 2280, column 127
Description:	The right operand of '*' is a garbage value

Annotated Source Code

#include "ANALYSIS/DHistogramActions.h"

void DHistogramAction_PID::Initialize(JEventLoop* locEventLoop)

{

string locHistName, locHistTitle;

string locParticleName, locParticleName2, locParticleROOTName, locParticleROOTName2;

Particle_t locPID, locPID2;

Run_Update(locEventLoop);

//auto locDesiredPIDs = Get_Reaction()->Get_FinalPIDs(-1, false, false, d_AllCharges, false);

auto locDesiredPIDs = Get_Reaction()->Get_FinalPIDs(-1, false, true, d_AllCharges, false);

vector<const DMCThrown*> locMCThrowns;

locEventLoop->Get(locMCThrowns);

//CREATE THE HISTOGRAMS

//Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock

japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!

{

CreateAndChangeTo_ActionDirectory();

for(size_t loc_i = 0; loc_i < locDesiredPIDs.size(); ++loc_i)

{

locPID = locDesiredPIDs[loc_i];

locParticleName = ParticleType(locPID);

locParticleROOTName = ParticleName_ROOT(locPID);

CreateAndChangeTo_Directory(locParticleName, locParticleName);

if(ParticleCharge(locPID) == 0 && Is_FinalStateParticle(locPID))

{

//BCAL

CreateAndChangeTo_Directory("BCAL", "BCAL");

locHistName = "Beta";

locHistTitle = string("BCAL ") + locParticleROOTName + string(" Candidates;#beta");

dHistMap_Beta[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumBetaBins, dMinBeta, dMaxBeta);

locHistName = "Shower_Energy";

locHistTitle = string("BCAL ") + locParticleROOTName + string(" Candidates;Shower Energy (GeV)");

dHistMap_CalE[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, 200, 0, 1.);

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

//FCAL

CreateAndChangeTo_Directory("FCAL", "FCAL");

locHistName = "Beta";

locHistTitle = string("FCAL ") + locParticleROOTName + string(" Candidates;#beta");

dHistMap_Beta[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumBetaBins, dMinBeta, dMaxBeta);

locHistName = "Shower_Energy";

locHistTitle = string("FCAL ") + locParticleROOTName + string(" Candidates;Shower Energy (GeV)");

dHistMap_CalE[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, 200, 0, 1.);

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

//CCAL

CreateAndChangeTo_Directory("CCAL", "CCAL");

locHistName = "Beta";

locHistTitle = string("CCAL ") + locParticleROOTName + string(" Candidates;#beta");

dHistMap_Beta[locPID][SYS_CCAL] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumBetaBins, dMinBeta, dMaxBeta);

locHistName = "Shower_Energy";

locHistTitle = string("CCAL ") + locParticleROOTName + string(" Candidates;Shower Energy (GeV)");

dHistMap_CalE[locPID][SYS_CCAL] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, 200, 0, 1.);

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

}

//q = 0

if(locPID == Gamma)

{

//BCAL

CreateAndChangeTo_Directory("BCAL", "BCAL");

locHistName = "BetaVsP";

locHistTitle = string("BCAL ") + locParticleROOTName + string(" Candidates;Shower Energy (GeV);#beta");

dHistMap_BetaVsP[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DBetaBins, dMinBeta, dMaxBeta);

locHistName = "DeltaTVsShowerE_Photon";

locHistTitle = string("BCAL ") + locParticleROOTName + string(" Candidates;Shower Energy (GeV);#Deltat_{BCAL - RF}");

dHistMap_DeltaTVsP[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);

locHistName = "TimePullVsShowerE_Photon";

locHistTitle = string("BCAL ") + locParticleROOTName + string(";Shower Energy (GeV);#Deltat/#sigma_{#Deltat}");

dHistMap_TimePullVsP[Gamma][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);

locHistName = "TimeFOMVsShowerE_Photon";

locHistTitle = string("BCAL ") + locParticleROOTName + string(";Shower Energy (GeV);Timing PID Confidence Level");

dHistMap_TimeFOMVsP[Gamma][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

//FCAL

CreateAndChangeTo_Directory("FCAL", "FCAL");

locHistName = "BetaVsP";

locHistTitle = string("FCAL ") + locParticleROOTName + string(" Candidates;Shower Energy (GeV);#beta");

dHistMap_BetaVsP[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);

locHistName = "DeltaTVsShowerE_Photon";

locHistTitle = string("FCAL ") + locParticleROOTName + string(";p (GeV/c);#Deltat_{FCAL - RF}");

dHistMap_DeltaTVsP[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);

locHistName = "TimePullVsShowerE_Photon";

locHistTitle = string("FCAL ") + locParticleROOTName + string(";Shower Energy (GeV);#Deltat/#sigma_{#Deltat}");

100

dHistMap_TimePullVsP[Gamma][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);

101

102

locHistName = "TimeFOMVsShowerE_Photon";

103

locHistTitle = string("FCAL ") + locParticleROOTName + string(";Shower Energy (GeV);Timing PID Confidence Level");

104

dHistMap_TimeFOMVsP[Gamma][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);

105

106

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

107

108

//CCAL

109

CreateAndChangeTo_Directory("CCAL", "CCAL");

110

111

locHistName = "BetaVsP";

112

locHistTitle = string("CCAL ") + locParticleROOTName + string(" Candidates;Shower Energy (GeV);#beta");

113

dHistMap_BetaVsP[locPID][SYS_CCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);

114

115

locHistName = "DeltaTVsShowerE_Photon";

116

locHistTitle = string("CCAL ") + locParticleROOTName + string(";p (GeV/c);#Deltat_{CCAL - RF}");

117

dHistMap_DeltaTVsP[locPID][SYS_CCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);

118

119

locHistName = "TimePullVsShowerE_Photon";

120

locHistTitle = string("CCAL ") + locParticleROOTName + string(";Shower Energy (GeV);#Deltat/#sigma_{#Deltat}");

121

dHistMap_TimePullVsP[Gamma][SYS_CCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);

122

123

locHistName = "TimeFOMVsShowerE_Photon";

124

locHistTitle = string("CCAL ") + locParticleROOTName + string(";Shower Energy (GeV);Timing PID Confidence Level");

125

dHistMap_TimeFOMVsP[Gamma][SYS_CCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);

126

127

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

128

}

129

else if(ParticleCharge(locPID) != 0)

130

{

131

//SC

132

CreateAndChangeTo_Directory("SC", "SC");

133

134

locHistName = "dEdXVsP";

135

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);SC dE/dX (MeV/cm)");

136

dHistMap_dEdXVsP[locPID][SYS_START] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);

137

138

locHistName = "DeltadEdXVsP";

139

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);SC #Delta(dE/dX) (MeV/cm)");

140

dHistMap_DeltadEdXVsP[locPID][SYS_START] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);

141

142

locHistName = "DeltaBetaVsP";

143

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);SC #Delta#beta");

144

dHistMap_DeltaBetaVsP[locPID][SYS_START] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);

145

146

locHistName = "BetaVsP";

147

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);SC #beta");

148

dHistMap_BetaVsP[locPID][SYS_START] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);

149

150

locHistName = "DeltaTVsP";

151

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{SC - RF}");

152

dHistMap_DeltaTVsP[locPID][SYS_START] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);

153

154

locHistName = string("TimePullVsP_") + locParticleName;

155

locHistTitle = string("SC ") + locParticleROOTName + string(";p (GeV/c);#Deltat/#sigma_{#Deltat}");

156

dHistMap_TimePullVsP[locPID][SYS_START] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);

157

158

locHistName = string("TimeFOMVsP_") + locParticleName;

159

locHistTitle = string("SC") + locParticleROOTName + string(";p (GeV/c);Timing PID Confidence Level");

160

dHistMap_TimeFOMVsP[locPID][SYS_START] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);

161

162

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

163

164

//TOF

165

CreateAndChangeTo_Directory("TOF", "TOF");

166

167

locHistName = "dEdXVsP";

168

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);TOF dE/dX (MeV/cm)");

169

dHistMap_dEdXVsP[locPID][SYS_TOF] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);

170

171

locHistName = "DeltadEdXVsP";

172

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);TOF #Delta(dE/dX) (MeV/cm)");

173

dHistMap_DeltadEdXVsP[locPID][SYS_TOF] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);

174

175

locHistName = "DeltaBetaVsP";

176

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);TOF #Delta#beta");

177

dHistMap_DeltaBetaVsP[locPID][SYS_TOF] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);

178

179

locHistName = "BetaVsP";

180

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);TOF #beta");

181

dHistMap_BetaVsP[locPID][SYS_TOF] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);

182

183

locHistName = "DeltaTVsP";

184

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{TOF - RF}");

185

dHistMap_DeltaTVsP[locPID][SYS_TOF] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);

186

187

locHistName = string("TimePullVsP_") + locParticleName;

188

locHistTitle = string("TOF ") + locParticleROOTName + string(";p (GeV/c);#Deltat/#sigma_{#Deltat}");

189

dHistMap_TimePullVsP[locPID][SYS_TOF] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);

190

191

locHistName = string("TimeFOMVsP_") + locParticleName;

192

locHistTitle = string("TOF ") + locParticleROOTName + string(";p (GeV/c);Timing PID Confidence Level");

193

dHistMap_TimeFOMVsP[locPID][SYS_TOF] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);

194

195

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

196

197

//BCAL

198

CreateAndChangeTo_Directory("BCAL", "BCAL");

199

200

locHistName = "BetaVsP";

201

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);BCAL #beta");

202

dHistMap_BetaVsP[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DBetaBins, dMinBeta, dMaxBeta);

203

204

locHistName = "DeltaBetaVsP";

205

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);BCAL #Delta#beta");

206

dHistMap_DeltaBetaVsP[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);

207

208

locHistName = "DeltaTVsP";

209

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{BCAL - RF}");

210

dHistMap_DeltaTVsP[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);

211

212

locHistName = string("TimePullVsP_") + locParticleName;

213

locHistTitle = string("BCAL ") + locParticleROOTName + string(";p (GeV/c);#Deltat/#sigma_{#Deltat}");

214

dHistMap_TimePullVsP[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);

215

216

locHistName = string("TimeFOMVsP_") + locParticleName;

217

locHistTitle = string("BCAL ") + locParticleROOTName + string(";p (GeV/c);Timing PID Confidence Level");

218

dHistMap_TimeFOMVsP[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);

219

220

locHistName = "EOverPVsP";

221

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);BCAL E_{Shower}/p_{Track} (c);");

222

dHistMap_EOverPVsP[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxBCALP, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);

223

224

locHistName = "EOverPVsTheta";

225

locHistTitle = locParticleROOTName + string(" Candidates;#theta#circ;BCAL E_{Shower}/p_{Track} (c);");

226

dHistMap_EOverPVsTheta[locPID][SYS_BCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DBCALThetaBins, dMinBCALTheta, dMaxBCALTheta, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);

227

228

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

229

230

//FCAL

231

CreateAndChangeTo_Directory("FCAL", "FCAL");

232

233

locHistName = "BetaVsP";

234

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);FCAL #beta");

235

dHistMap_BetaVsP[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);

236

237

locHistName = "DeltaBetaVsP";

238

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);FCAL #Delta#beta");

239

dHistMap_DeltaBetaVsP[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);

240

241

locHistName = "DeltaTVsP";

242

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{FCAL - RF}");

243

dHistMap_DeltaTVsP[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);

244

245

locHistName = string("TimePullVsP_") + locParticleName;

246

locHistTitle = string("FCAL ") + locParticleROOTName + string(";p (GeV/c);#Deltat/#sigma_{#Deltat}");

247

dHistMap_TimePullVsP[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);

248

249

locHistName = string("TimeFOMVsP_") + locParticleName;

250

locHistTitle = string("FCAL ") + locParticleROOTName + string(";p (GeV/c);Timing PID Confidence Level");

251

dHistMap_TimeFOMVsP[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);

252

253

locHistName = "EOverPVsP";

254

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);FCAL E_{Shower}/p_{Track} (c);");

255

dHistMap_EOverPVsP[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);

256

257

locHistName = "EOverPVsTheta";

258

locHistTitle = locParticleROOTName + string(" Candidates;#theta#circ;FCAL E_{Shower}/p_{Track} (c);");

259

dHistMap_EOverPVsTheta[locPID][SYS_FCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DFCALThetaBins, dMinFCALTheta, dMaxFCALTheta, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);

260

261

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

262

263

//CCAL

264

CreateAndChangeTo_Directory("CCAL", "CCAL");

265

266

locHistName = "BetaVsP";

267

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CCAL #beta");

268

dHistMap_BetaVsP[locPID][SYS_CCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DBetaBins, dMinBeta, dMaxBeta);

269

270

locHistName = "DeltaBetaVsP";

271

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CCAL #Delta#beta");

272

dHistMap_DeltaBetaVsP[locPID][SYS_CCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);

273

274

locHistName = "DeltaTVsP";

275

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);#Deltat_{CCAL - RF}");

276

dHistMap_DeltaTVsP[locPID][SYS_CCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltaTBins, dMinDeltaT, dMaxDeltaT);

277

278

locHistName = string("TimePullVsP_") + locParticleName;

279

locHistTitle = string("CCAL ") + locParticleROOTName + string(";p (GeV/c);#Deltat/#sigma_{#Deltat}");

280

dHistMap_TimePullVsP[locPID][SYS_CCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);

281

282

locHistName = string("TimeFOMVsP_") + locParticleName;

283

locHistTitle = string("CCAL ") + locParticleROOTName + string(";p (GeV/c);Timing PID Confidence Level");

284

dHistMap_TimeFOMVsP[locPID][SYS_CCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);

285

286

locHistName = "EOverPVsP";

287

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CCAL E_{Shower}/p_{Track} (c);");

288

dHistMap_EOverPVsP[locPID][SYS_CCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);

289

290

locHistName = "EOverPVsTheta";

291

locHistTitle = locParticleROOTName + string(" Candidates;#theta#circ;CCAL E_{Shower}/p_{Track} (c);");

292

dHistMap_EOverPVsTheta[locPID][SYS_CCAL] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DCCALThetaBins, dMinCCALTheta, dMaxCCALTheta, dNum2DEOverPBins, dMinEOverP, dMaxEOverP);

293

294

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

295

296

//CDC

297

CreateAndChangeTo_Directory("CDC", "CDC");

298

299

locHistName = "dEdXVsP_Int";

300

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CDC dE/dx [Integral-based] (keV/cm)");

301

dHistMap_dEdXVsP[locPID][SYS_CDC] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);

302

303

locHistName = "dEdXVsP_Amp";

304

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CDC dE/dx [Amplitude-based] (keV/cm)");

305

dHistMap_dEdXVsP[locPID][SYS_CDC_AMP] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);

306

307

locHistName = "DeltadEdXVsP_Int";

308

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CDC #Delta(dE/dX) [Integral-based] (keV/cm)");

309

dHistMap_DeltadEdXVsP[locPID][SYS_CDC] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);

310

311

locHistName = "DeltadEdXVsP_Amp";

312

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);CDC #Delta(dE/dX) [Amplitude-based] (keV/cm)");

313

dHistMap_DeltadEdXVsP[locPID][SYS_CDC_AMP] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);

314

315

locHistName = string("dEdXPullVsP_Int_") + locParticleName;

316

locHistTitle = locParticleROOTName + string(";p (GeV/c);CDC #Delta(dE/dX)/#sigma_{#Delta(dE/dX)} [Integral-based]");

317

dHistMap_dEdXPullVsP[locPID][SYS_CDC] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);

318

319

locHistName = string("dEdXPullVsP_Amp_") + locParticleName;

320

locHistTitle = locParticleROOTName + string(";p (GeV/c);CDC #Delta(dE/dX)/#sigma_{#Delta(dE/dX)} [Amplitude-based]");

321

dHistMap_dEdXPullVsP[locPID][SYS_CDC_AMP] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);

322

323

locHistName = string("dEdXFOMVsP_Int_") + locParticleName;

324

locHistTitle = locParticleROOTName + string(";p (GeV/c);CDC dE/dx [Integral-based] PID Confidence Level");

325

dHistMap_dEdXFOMVsP[locPID][SYS_CDC] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);

326

327

locHistName = string("dEdXFOMVsP_Amp_") + locParticleName;

328

locHistTitle = locParticleROOTName + string(";p (GeV/c);CDC dE/dx [Amplitude-based] PID Confidence Level");

329

dHistMap_dEdXFOMVsP[locPID][SYS_CDC_AMP] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);

330

331

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

332

333

//FDC

334

CreateAndChangeTo_Directory("FDC", "FDC");

335

336

locHistName = "dEdXVsP";

337

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);FDC dE/dx (keV/cm)");

338

dHistMap_dEdXVsP[locPID][SYS_FDC] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DdEdxBins, dMindEdX, dMaxdEdX);

339

340

locHistName = "DeltadEdXVsP";

341

locHistTitle = locParticleROOTName + string(" Candidates;p (GeV/c);FDC #Delta(dE/dX) (keV/cm)");

342

dHistMap_DeltadEdXVsP[locPID][SYS_FDC] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DDeltadEdxBins, dMinDeltadEdx, dMaxDeltadEdx);

343

344

locHistName = string("dEdXPullVsP_") + locParticleName;

345

locHistTitle = locParticleROOTName + string(";p (GeV/c);FDC #Delta(dE/dX)/#sigma_{#Delta(dE/dX)}");

346

dHistMap_dEdXPullVsP[locPID][SYS_FDC] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DPullBins, dMinPull, dMaxPull);

347

348

locHistName = string("dEdXFOMVsP_") + locParticleName;

349

locHistTitle = locParticleROOTName + string(";p (GeV/c);FDC dE/dx PID Confidence Level");

350

dHistMap_dEdXFOMVsP[locPID][SYS_FDC] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNum2DFOMBins, 0.0, 1.0);

351

352

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

353

354

// DIRC

355

CreateAndChangeTo_Directory("DIRC", "DIRC");

356

357

locHistName = string("NumPhotons_") + locParticleName;

358

locHistTitle = locParticleROOTName + string("; DIRC NumPhotons");

359

dHistMap_NumPhotons_DIRC[locPID] = new TH1I(locHistName.c_str(), locHistTitle.c_str(), dDIRCNumPhotonsBins, dDIRCMinNumPhotons, dDIRCMaxNumPhotons);

360

361

locHistName = string("ThetaCVsP_") + locParticleName;

362

locHistTitle = locParticleROOTName + string("; Momentum (GeV); DIRC #theta_{C}");

363

dHistMap_ThetaCVsP_DIRC[locPID] = new TH2I(locHistName.c_str(), locHistTitle.c_str(), dNum2DPBins, dMinP, dMaxP, dDIRCThetaCBins, dDIRCMinThetaC, dDIRCMaxThetaC);

364

365

locHistName = string("Ldiff_kpiVsP_") + locParticleName;

366

locHistTitle = locParticleROOTName + string("; Momentum (GeV); DIRC L_{K}-L_{#pi}");

367

dHistMap_Ldiff_kpiVsP_DIRC[locPID] = new TH2I(locHistName.c_str(), locHistTitle.c_str(), dNum2DPBins, dMinP, dMaxP, dDIRCLikelihoodBins, -1*dDIRCMaxLikelihood, dDIRCMaxLikelihood);

368

369

locHistName = string("Ldiff_pkVsP_") + locParticleName;

370

locHistTitle = locParticleROOTName + string("; Momentum (GeV); DIRC L_{p}-L_{K}");

371

dHistMap_Ldiff_pkVsP_DIRC[locPID] = new TH2I(locHistName.c_str(), locHistTitle.c_str(), dNum2DPBins, dMinP, dMaxP, dDIRCLikelihoodBins, -1*dDIRCMaxLikelihood, dDIRCMaxLikelihood);

372

373

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

374

} //end of charged

375

376

//one per thrown pid:

377

if(!locMCThrowns.empty())

378

{

379

CreateAndChangeTo_Directory("Throwns", "Throwns");

380

for(size_t loc_j = 0; loc_j < dThrownPIDs.size(); ++loc_j)

381

{

382

locPID2 = dThrownPIDs[loc_j];

383

if((ParticleCharge(locPID2) != ParticleCharge(locPID)) && (locPID2 != Unknown))

384

continue;

385

locParticleName2 = ParticleType(locPID2);

386

locParticleROOTName2 = ParticleName_ROOT(locPID2);

387

388

//Confidence Level for Thrown PID

389

pair<Particle_t, Particle_t> locPIDPair(locPID, locPID2);

390

locHistName = string("PIDConfidenceLevel_ForThrownPID_") + locParticleName2;

391

locHistTitle = locParticleROOTName + string(" PID, Thrown PID = ") + locParticleROOTName2 + string(";PID Confidence Level");

392

dHistMap_PIDFOMForTruePID[locPIDPair] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumFOMBins, 0.0, 1.0);

393

}

394

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

395

}

396

397

//no FOM for massive neutrals

398

if((ParticleCharge(locPID) != 0 && Is_FinalStateParticle(locPID)) || (locPID == Gamma))

399

{

400

// Overall Confidence Level

401

locHistName = "PIDConfidenceLevel";

402

locHistTitle = locParticleROOTName + string(" PID;PID Confidence Level");

403

dHistMap_PIDFOM[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumFOMBins, 0.0, 1.0);

404

405

// P Vs Theta, PID FOM = NaN

406

locHistName = "PVsTheta_NaNPIDFOM";

407

locHistTitle = locParticleROOTName + string(", PID FOM = NaN;#theta#circ;p (GeV/c)");

408

dHistMap_PVsTheta_NaNPIDFOM[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);

409

}

410

411

// look at decaying particles

412

if(!Is_FinalStateParticle(locPID))

413

{

414

// Flight Distance

415

locHistName = "FlightDistance";

416

locHistTitle = locParticleROOTName + string(" Flight Distance;Flight Distance (cm)");

417

dHistMap_FlightDistance[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumFlightDistanceBins, 0.0, dMaxFlightDistance);

418

419

// Flight Distance Vs P

420

locHistName = "FlightDistanceVsP";

421

locHistTitle = locParticleROOTName + string(" Flight Distance;Flight Distance (cm);p (GeV/c)");

422

dHistMap_FlightDistanceVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumFlightDistanceBins, 0.0, dMaxFlightDistance, dNum2DPBins, dMinP, dMaxP);

423

424

// Flight Distance Vs Theta

425

locHistName = "FlightDistanceVsTheta";

426

locHistTitle = locParticleROOTName + string(" Flight Distance;Flight Distance (cm);#theta#circ");

427

dHistMap_FlightDistanceVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumFlightDistanceBins, 0.0, dMaxFlightDistance, dNum2DThetaBins, dMinTheta, dMaxTheta);

428

429

430

// Flight Distance

431

locHistName = "FlightSignificance";

432

locHistTitle = locParticleROOTName + string(" Flight Significance;Flight Significance (#sigma)");

433

dHistMap_FlightSignificance[locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumFlightSignificanceBins, 0.0, dMaxFlightSignificance);

434

435

// Flight Distance Vs P

436

locHistName = "FlightSignificanceVsP";

437

locHistTitle = locParticleROOTName + string(" Flight Significance;Flight Significance (#sigma);p (GeV/c)");

438

dHistMap_FlightSignificanceVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumFlightSignificanceBins, 0.0, dMaxFlightSignificance, dNum2DPBins, dMinP, dMaxP);

439

440

// Flight Distance Vs Theta

441

locHistName = "FlightSignificanceVsTheta";

442

locHistTitle = locParticleROOTName + string(" Flight Significance;Flight Significance (#sigma);#theta#circ");

443

dHistMap_FlightSignificanceVsTheta[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumFlightSignificanceBins, 0.0, dMaxFlightSignificance, dNum2DThetaBins, dMinTheta, dMaxTheta);

444

445

}

446

447

448

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

449

} //end of PID loop

450

451

//Return to the base directory

452

ChangeTo_BaseDirectory();

453

}

454

japp->RootUnLock(); //RELEASE ROOT LOCK!!

455

}

456

457

void DHistogramAction_PID::Run_Update(JEventLoop* locEventLoop)

458

{

459

vector<const DParticleID*> locParticleIDs;

460

locEventLoop->Get(locParticleIDs);

461

462

vector<const DAnalysisUtilities*> locAnalysisUtilitiesVector;

463

locEventLoop->Get(locAnalysisUtilitiesVector);

464

465

dParticleID = locParticleIDs[0];

466

dAnalysisUtilities = locAnalysisUtilitiesVector[0];

467

}

468

469

bool DHistogramAction_PID::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)

470

{

471

vector<const DMCThrownMatching*> locMCThrownMatchingVector;

472

locEventLoop->Get(locMCThrownMatchingVector);

473

const DMCThrownMatching* locMCThrownMatching = locMCThrownMatchingVector.empty() ? NULL__null : locMCThrownMatchingVector[0];

474

475

const DEventRFBunch* locEventRFBunch = locParticleCombo->Get_EventRFBunch();

476

477

vector<const DReactionVertexInfo*> locVertexInfos;

478

locEventLoop->Get(locVertexInfos);

479

480

// figure out what the right DReactionVertexInfo is

481

const DReactionVertexInfo *locReactionVertexInfo = nullptr;

482

for(auto& locVertexInfo : locVertexInfos) {

483

auto locVertexReactions = locVertexInfo->Get_Reactions();

484

if(find(locVertexReactions.begin(), locVertexReactions.end(), Get_Reaction()) != locVertexReactions.end()) {

485

locReactionVertexInfo = locVertexInfo;

486

break;

487

}

488

}

489

490

for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i)

491

{

492

const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(loc_i);

493

auto locParticles = Get_UseKinFitResultsFlag() ? locParticleComboStep->Get_FinalParticles(Get_Reaction()->Get_ReactionStep(loc_i), false, false, d_AllCharges) : locParticleComboStep->Get_FinalParticles_Measured(Get_Reaction()->Get_ReactionStep(loc_i), d_AllCharges);

494

for(size_t loc_j = 0; loc_j < locParticles.size(); ++loc_j)

495

{

496

//check if will be duplicate

497

pair<Particle_t, const JObject*> locParticleInfo(locParticles[loc_j]->PID(), Get_FinalParticle_SourceObject(locParticles[loc_j]));

498

pair<const DEventRFBunch*, pair<Particle_t, const JObject*> > locHistInfo(locEventRFBunch, locParticleInfo);

499

if(dPreviouslyHistogrammedParticles.find(locHistInfo) != dPreviouslyHistogrammedParticles.end())

500

continue; //previously histogrammed

501

502

dPreviouslyHistogrammedParticles.insert(locHistInfo);

503

if(Is_FinalStateParticle(locParticles[loc_j]->PID())) {

504

if(ParticleCharge(locParticles[loc_j]->PID()) != 0) //charged

505

Fill_ChargedHists(static_cast<const DChargedTrackHypothesis*>(locParticles[loc_j]), locMCThrownMatching, locEventRFBunch);

506

else //neutral

507

Fill_NeutralHists(static_cast<const DNeutralParticleHypothesis*>(locParticles[loc_j]), locMCThrownMatching, locEventRFBunch);

508

}

509

}

510

511

// fill info on decaying particles, which is on the particle combo step level

512

if(((locParticleComboStep->Get_InitialParticle()!=nullptr) && !Is_FinalStateParticle(locParticleComboStep->Get_InitialParticle()->PID()))) { // decaying particles

513

// FOR NOW: we only calculate these displaced vertex quantities if a kinematic fit

514

// was performed where the vertex is constrained

515

// TODO: Cleverly handle the other cases

516

517

// pull out information related to the kinematic fit

518

if( !Get_UseKinFitResultsFlag() ) continue;

519

if( locReactionVertexInfo == nullptr ) continue;

520

521

// extract the info about the vertex, to make sure that the information that we

522

// need to calculate displaced vertices is there

523

auto locStepVertexInfo = locReactionVertexInfo->Get_StepVertexInfo(loc_i);

524

525

Fill_DecayingHists(locParticleComboStep->Get_InitialParticle(), locParticleComboStep->Get_InitialKinFitParticle(), locStepVertexInfo);

526

}

527

528

}

529

530

return true;

531

}

532

533

void DHistogramAction_PID::Fill_DecayingHists(const DKinematicData *locInitialParticle, std::shared_ptr<const DKinFitParticle> locKinFitParticle, const DReactionStepVertexInfo *locStepVertexInfo)

534

{

535

DKinFitType locKinFitType = Get_Reaction()->Get_KinFitType();

536

537

//Particle_t locPID = static_cast<Particle_t>(locKinFitParticle->Get_PID());

538

//double locP = locKinFitParticle->Get_Momentum().Mag();

539

//double locTheta = locKinFitParticle->Get_Momentum().Theta()*180.0/TMath::Pi();

540

Particle_t locPID = locInitialParticle->PID();

541

double locP = locInitialParticle->momentum().Mag();

542

double locTheta = locInitialParticle->momentum().Theta()*180.0/TMath::Pi();

543

544

// comment

545

auto locParentVertexInfo = locStepVertexInfo->Get_ParentVertexInfo();

546

auto locVertexKinFitFlag = ((locKinFitType != d_P4Fit) && (locKinFitType != d_NoFit));

547

548

if(IsDetachedVertex(locPID) && locVertexKinFitFlag && (locParentVertexInfo != nullptr) && locStepVertexInfo->Get_FittableVertexFlag() && locParentVertexInfo->Get_FittableVertexFlag())

549

{

550

auto locPathLength = locKinFitParticle->Get_PathLength();

551

auto locPathLengthSigma = locKinFitParticle->Get_PathLengthUncertainty();

552

double locPathLengthSignificance = locPathLength / locPathLengthSigma;

553

554

dHistMap_FlightDistance[locPID]->Fill(locPathLength);

555

dHistMap_FlightDistanceVsP[locPID]->Fill(locPathLength, locP);

556

dHistMap_FlightDistanceVsTheta[locPID]->Fill(locPathLength, locTheta);

557

558

dHistMap_FlightSignificance[locPID]->Fill(locPathLengthSignificance);

559

dHistMap_FlightSignificanceVsP[locPID]->Fill(locPathLengthSignificance, locP);

560

dHistMap_FlightSignificanceVsTheta[locPID]->Fill(locPathLengthSignificance, locTheta);

561

562

}

563

564

}

565

566

567

void DHistogramAction_PID::Fill_ChargedHists(const DChargedTrackHypothesis* locChargedTrackHypothesis, const DMCThrownMatching* locMCThrownMatching, const DEventRFBunch* locEventRFBunch)

568

{

569

Particle_t locPID = locChargedTrackHypothesis->PID();

570

double locBeta_Timing = locChargedTrackHypothesis->measuredBeta();

571

double locDeltaBeta = locBeta_Timing - locChargedTrackHypothesis->lorentzMomentum().Beta();

572

double locDeltaT = (locChargedTrackHypothesis->Get_TimeAtPOCAToVertex() - locChargedTrackHypothesis->t0());

573

double locFOM_Timing = (locChargedTrackHypothesis->Get_NDF_Timing() > 0) ? TMath::Prob(locChargedTrackHypothesis->Get_ChiSq_Timing(), locChargedTrackHypothesis->Get_NDF_Timing()) : numeric_limits<double>::quiet_NaN();

574

575

double locP = locChargedTrackHypothesis->momentum().Mag();

576

double locTheta = locChargedTrackHypothesis->momentum().Theta()*180.0/TMath::Pi();

577

double locMatchFOM = 0.0;

578

const DMCThrown* locMCThrown = (locMCThrownMatching != NULL__null) ? locMCThrownMatching->Get_MatchingMCThrown(locChargedTrackHypothesis, locMatchFOM) : NULL__null;

579

580

auto locBCALShowerMatchParams = locChargedTrackHypothesis->Get_BCALShowerMatchParams();

581

auto locFCALShowerMatchParams = locChargedTrackHypothesis->Get_FCALShowerMatchParams();

582

auto locTOFHitMatchParams = locChargedTrackHypothesis->Get_TOFHitMatchParams();

583

auto locSCHitMatchParams = locChargedTrackHypothesis->Get_SCHitMatchParams();

584

auto locDIRCMatchParams = locChargedTrackHypothesis->Get_DIRCMatchParams();

585

586

auto locTrackTimeBased = locChargedTrackHypothesis->Get_TrackTimeBased();

587

588

double locTimePull = 0.0;

589

unsigned int locTimeNDF = 0;

590

dParticleID->Calc_TimingChiSq(locChargedTrackHypothesis, locTimeNDF, locTimePull);

591

DetectorSystem_t locTimeDetector = locChargedTrackHypothesis->t1_detector();

592

593

//FILL HISTOGRAMS

594

//Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock

595

//Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex

596

Lock_Action();

597

{

598

dHistMap_PIDFOM[locPID]->Fill(locChargedTrackHypothesis->Get_FOM());

599

600

//SC dE/dx

601

if(locSCHitMatchParams != NULL__null)

602

{

603

dHistMap_dEdXVsP[locPID][SYS_START]->Fill(locP, locSCHitMatchParams->dEdx*1.0E3);

604

double locdx = locSCHitMatchParams->dHitEnergy/locSCHitMatchParams->dEdx;

605

double locProbabledEdx = 0.0, locSigmadEdx = 0.0;

606

dParticleID->GetScintMPdEandSigma(locP, locChargedTrackHypothesis->mass(), locdx, locProbabledEdx, locSigmadEdx);

607

dHistMap_DeltadEdXVsP[locPID][SYS_START]->Fill(locP, (locSCHitMatchParams->dEdx - locProbabledEdx)*1.0E3);

608

}

609

610

//TOF dE/dx

611

if(locTOFHitMatchParams != NULL__null)

612

{

613

dHistMap_dEdXVsP[locPID][SYS_TOF]->Fill(locP, locTOFHitMatchParams->dEdx*1.0E3);

614

double locdx = locTOFHitMatchParams->dHitEnergy/locTOFHitMatchParams->dEdx;

615

double locProbabledEdx = 0.0, locSigmadEdx = 0.0;

616

dParticleID->GetScintMPdEandSigma(locP, locChargedTrackHypothesis->mass(), locdx, locProbabledEdx, locSigmadEdx);

617

dHistMap_DeltadEdXVsP[locPID][SYS_TOF]->Fill(locP, (locTOFHitMatchParams->dEdx - locProbabledEdx)*1.0E3);

618

}

619

620

//BCAL E/p

621

if(locBCALShowerMatchParams != NULL__null)

622

{

623

const DBCALShower* locBCALShower = locBCALShowerMatchParams->dBCALShower;

624

double locEOverP = locBCALShower->E/locP;

625

dHistMap_EOverPVsP[locPID][SYS_BCAL]->Fill(locP, locEOverP);

626

dHistMap_EOverPVsTheta[locPID][SYS_BCAL]->Fill(locTheta, locEOverP);

627

}

628

629

//FCAL E/p

630

if(locFCALShowerMatchParams != NULL__null)

631

{

632

const DFCALShower* locFCALShower = locFCALShowerMatchParams->dFCALShower;

633

double locEOverP = locFCALShower->getEnergy()/locP;

634

dHistMap_EOverPVsP[locPID][SYS_FCAL]->Fill(locP, locEOverP);

635

dHistMap_EOverPVsTheta[locPID][SYS_FCAL]->Fill(locTheta, locEOverP);

636

}

637

638

//CCAL E/p

639

if(locCCALShowerMatchParams != NULL)

640

{

641

const DFCALShower* locCCALShower = locCCALShowerMatchParams->dFCALShower;

642

double locEOverP = locCCALShower->getEnergy()/locP;

643

dHistMap_EOverPVsP[locPID][SYS_CCAL]->Fill(locP, locEOverP);

644

dHistMap_EOverPVsTheta[locPID][SYS_CCAL]->Fill(locTheta, locEOverP);

645

}

646

647

//DIRC

648

if(locDIRCMatchParams != NULL__null) {

649

650

int locNumPhotons_DIRC = locDIRCMatchParams->dNPhotons;

651

double locThetaC_DIRC = locDIRCMatchParams->dThetaC;

652

dHistMap_NumPhotons_DIRC[locPID]->Fill(locNumPhotons_DIRC);

653

dHistMap_ThetaCVsP_DIRC[locPID]->Fill(locP, locThetaC_DIRC);

654

double locLpi_DIRC = locDIRCMatchParams->dLikelihoodPion;

655

double locLk_DIRC = locDIRCMatchParams->dLikelihoodKaon;

656

double locLp_DIRC = locDIRCMatchParams->dLikelihoodProton;

657

dHistMap_Ldiff_kpiVsP_DIRC[locPID]->Fill(locP, locLk_DIRC-locLpi_DIRC);

658

dHistMap_Ldiff_pkVsP_DIRC[locPID]->Fill(locP, locLp_DIRC-locLk_DIRC);

659

}

660

661

//Timing

662

if(dHistMap_BetaVsP[locPID].find(locTimeDetector) != dHistMap_BetaVsP[locPID].end())

663

{

664

dHistMap_BetaVsP[locPID][locTimeDetector]->Fill(locP, locBeta_Timing);

665

dHistMap_DeltaBetaVsP[locPID][locTimeDetector]->Fill(locP, locDeltaBeta);

666

dHistMap_DeltaTVsP[locPID][locTimeDetector]->Fill(locP, locDeltaT);

667

dHistMap_TimePullVsP[locPID][locTimeDetector]->Fill(locP, locTimePull);

668

dHistMap_TimeFOMVsP[locPID][locTimeDetector]->Fill(locP, locFOM_Timing);

669

}

670

671

//CDC dE/dx

672

if(locTrackTimeBased->dNumHitsUsedFordEdx_CDC > 0)

673

{

674

dHistMap_dEdXVsP[locPID][SYS_CDC]->Fill(locP, locChargedTrackHypothesis->Get_dEdx_CDC_int()*1.0E6);

675

dHistMap_dEdXVsP[locPID][SYS_CDC_AMP]->Fill(locP, locChargedTrackHypothesis->Get_dEdx_CDC_amp()*1.0E6);

676

677

double locProbabledEdx = dParticleID->GetMostProbabledEdx_DC(locP, locChargedTrackHypothesis->mass(), locTrackTimeBased->ddx_CDC, true);

678

double locDeltadEdx = locChargedTrackHypothesis->Get_dEdx_CDC_int() - locProbabledEdx;

679

dHistMap_DeltadEdXVsP[locPID][SYS_CDC]->Fill(locP, 1.0E6*locDeltadEdx);

680

double locDeltadEdx_amp = locChargedTrackHypothesis->Get_dEdx_CDC_amp() - locProbabledEdx;

681

dHistMap_DeltadEdXVsP[locPID][SYS_CDC_AMP]->Fill(locP, 1.0E6*locDeltadEdx_amp);

682

683

double locMeandx = locTrackTimeBased->ddx_CDC/locTrackTimeBased->dNumHitsUsedFordEdx_CDC;

684

double locSigmadEdx = dParticleID->GetdEdxSigma_DC(locTrackTimeBased->dNumHitsUsedFordEdx_CDC, locP, locChargedTrackHypothesis->mass(), locMeandx, true);

685

double locdEdXPull = locDeltadEdx/locSigmadEdx;

686

dHistMap_dEdXPullVsP[locPID][SYS_CDC]->Fill(locP, locDeltadEdx/locSigmadEdx);

687

double locdEdXPull_amp = locDeltadEdx_amp/locSigmadEdx;

688

dHistMap_dEdXPullVsP[locPID][SYS_CDC_AMP]->Fill(locP, locDeltadEdx_amp/locSigmadEdx);

689

690

double locdEdXChiSq = locdEdXPull*locdEdXPull;

691

double locdEdXFOM = TMath::Prob(locdEdXChiSq, locTrackTimeBased->dNumHitsUsedFordEdx_CDC);

692

dHistMap_dEdXFOMVsP[locPID][SYS_CDC]->Fill(locP, locdEdXFOM);

693

double locdEdXChiSq_amp = locdEdXPull_amp*locdEdXPull_amp;

694

double locdEdXFOM_amp = TMath::Prob(locdEdXChiSq_amp, locTrackTimeBased->dNumHitsUsedFordEdx_CDC);

695

dHistMap_dEdXFOMVsP[locPID][SYS_CDC_AMP]->Fill(locP, locdEdXFOM_amp);

696

}

697

698

//FDC dE/dx

699

if(locTrackTimeBased->dNumHitsUsedFordEdx_FDC > 0)

700

{

701

dHistMap_dEdXVsP[locPID][SYS_FDC]->Fill(locP, locTrackTimeBased->ddEdx_FDC*1.0E6);

702

double locProbabledEdx = dParticleID->GetMostProbabledEdx_DC(locP, locChargedTrackHypothesis->mass(), locTrackTimeBased->ddx_FDC, false);

703

double locDeltadEdx = locTrackTimeBased->ddEdx_FDC - locProbabledEdx;

704

dHistMap_DeltadEdXVsP[locPID][SYS_FDC]->Fill(locP, 1.0E6*locDeltadEdx);

705

double locMeandx = locTrackTimeBased->ddx_FDC/locTrackTimeBased->dNumHitsUsedFordEdx_FDC;

706

double locSigmadEdx = dParticleID->GetdEdxSigma_DC(locTrackTimeBased->dNumHitsUsedFordEdx_FDC, locP, locChargedTrackHypothesis->mass(), locMeandx, false);

707

double locdEdXPull = locDeltadEdx/locSigmadEdx;

708

dHistMap_dEdXPullVsP[locPID][SYS_FDC]->Fill(locP, locDeltadEdx/locSigmadEdx);

709

double locdEdXChiSq = locdEdXPull*locdEdXPull;

710

double locdEdXFOM = TMath::Prob(locdEdXChiSq, locTrackTimeBased->dNumHitsUsedFordEdx_FDC);

711

dHistMap_dEdXFOMVsP[locPID][SYS_FDC]->Fill(locP, locdEdXFOM);

712

}

713

714

pair<Particle_t, Particle_t> locPIDPair(locPID, Unknown); //default unless matched

715

if(locMCThrown != NULL__null) //else bogus track (not matched to any thrown tracks)

716

locPIDPair.second = (Particle_t)(locMCThrown->type); //matched

717

if(dHistMap_PIDFOMForTruePID.find(locPIDPair) != dHistMap_PIDFOMForTruePID.end()) //else hist not created or PID is weird

718

dHistMap_PIDFOMForTruePID[locPIDPair]->Fill(locChargedTrackHypothesis->Get_FOM());

719

720

if(locChargedTrackHypothesis->Get_NDF() == 0) //NaN

721

dHistMap_PVsTheta_NaNPIDFOM[locPID]->Fill(locTheta, locP);

722

}

723

Unlock_Action();

724

}

725

726

void DHistogramAction_PID::Fill_NeutralHists(const DNeutralParticleHypothesis* locNeutralParticleHypothesis, const DMCThrownMatching* locMCThrownMatching, const DEventRFBunch* locEventRFBunch)

727

{

728

Particle_t locPID = locNeutralParticleHypothesis->PID();

729

730

double locBeta_Timing = locNeutralParticleHypothesis->measuredBeta();

731

double locDeltaT = locNeutralParticleHypothesis->time() - locNeutralParticleHypothesis->t0();

732

733

double locP = locNeutralParticleHypothesis->momentum().Mag();

734

double locShowerE = locNeutralParticleHypothesis->Get_NeutralShower()->dEnergy;

735

double locMatchFOM = 0.0;

736

const DMCThrown* locMCThrown = (locMCThrownMatching != NULL__null) ? locMCThrownMatching->Get_MatchingMCThrown(locNeutralParticleHypothesis, locMatchFOM) : NULL__null;

737

738

double locTimePull = 0.0;

739

unsigned int locTimeNDF = 0;

740

dParticleID->Calc_TimingChiSq(locNeutralParticleHypothesis, locTimeNDF, locTimePull);

741

742

DetectorSystem_t locSystem = locNeutralParticleHypothesis->t1_detector();

743

744

//FILL HISTOGRAMS

745

//Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock

746

//Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex

747

Lock_Action();

748

{

749

//Beta (good for all PIDs)

750

dHistMap_Beta[locPID][locSystem]->Fill(locBeta_Timing);

751

dHistMap_CalE[locPID][locSystem]->Fill(locShowerE);

752

if(locPID != Gamma)

753

{

754

Unlock_Action();

755

return;

756

}

757

758

dHistMap_PIDFOM[locPID]->Fill(locNeutralParticleHypothesis->Get_FOM());

759

dHistMap_BetaVsP[locPID][locSystem]->Fill(locP, locBeta_Timing);

760

dHistMap_DeltaTVsP[locPID][locSystem]->Fill(locP, locDeltaT);

761

dHistMap_TimePullVsP[locPID][locSystem]->Fill(locP, locTimePull);

762

dHistMap_TimeFOMVsP[locPID][locSystem]->Fill(locP, locNeutralParticleHypothesis->Get_FOM());

763

764

pair<Particle_t, Particle_t> locPIDPair(locPID, Unknown); //default unless matched

765

if(locMCThrown != NULL__null) //else bogus track (not matched to any thrown tracks)

766

locPIDPair.second = (Particle_t)(locMCThrown->type); //matched

767

if(dHistMap_PIDFOMForTruePID.find(locPIDPair) != dHistMap_PIDFOMForTruePID.end()) //else hist not created or PID is weird

768

dHistMap_PIDFOMForTruePID[locPIDPair]->Fill(locNeutralParticleHypothesis->Get_FOM());

769

}

770

Unlock_Action();

771

}

772

773

void DHistogramAction_TrackVertexComparison::Initialize(JEventLoop* locEventLoop)

774

{

775

string locHistName, locHistTitle, locStepROOTName, locParticleName, locParticleROOTName;

776

Particle_t locPID, locHigherMassPID, locLowerMassPID;

777

string locHigherMassParticleName, locLowerMassParticleName, locHigherMassParticleROOTName, locLowerMassParticleROOTName;

778

779

Run_Update(locEventLoop);

780

781

size_t locNumSteps = Get_Reaction()->Get_NumReactionSteps();

782

dHistDeque_TrackZToCommon.resize(locNumSteps);

783

dHistDeque_TrackTToCommon.resize(locNumSteps);

784

dHistDeque_TrackDOCAToCommon.resize(locNumSteps);

785

dHistDeque_MaxTrackDeltaZ.resize(locNumSteps);

786

dHistDeque_MaxTrackDeltaT.resize(locNumSteps);

787

dHistDeque_MaxTrackDOCA.resize(locNumSteps);

788

dHistDeque_TrackDeltaTVsP.resize(locNumSteps);

789

790

//CREATE THE HISTOGRAMS

791

//Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock

792

japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!

793

{

794

CreateAndChangeTo_ActionDirectory();

795

for(size_t loc_i = 0; loc_i < locNumSteps; ++loc_i)

796

{

797

auto locDetectedChargedPIDs = Get_Reaction()->Get_FinalPIDs(loc_i, false, false, d_Charged, false);

798

if(locDetectedChargedPIDs.empty())

799

continue;

800

801

const DReactionStep* locReactionStep = Get_Reaction()->Get_ReactionStep(loc_i);

802

ostringstream locStepName;

803

locStepName << "Step" << loc_i << "__" << DAnalysis::Get_StepName(locReactionStep, true, false);

804

locStepROOTName = DAnalysis::Get_StepName(locReactionStep, true, true);

805

CreateAndChangeTo_Directory(locStepName.str(), locStepName.str());

806

807

// Max Track DeltaZ

808

locHistName = "MaxTrackDeltaZ";

809

locHistTitle = locStepROOTName + string(";Largest Track #DeltaVertex-Z (cm)");

810

dHistDeque_MaxTrackDeltaZ[loc_i] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaVertexZBins, 0.0, dMaxDeltaVertexZ);

811

812

// Max Track DeltaT

813

locHistName = "MaxTrackDeltaT";

814

locHistTitle = locStepROOTName + string(";Largest Track #DeltaVertex-T (ns)");

815

dHistDeque_MaxTrackDeltaT[loc_i] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaVertexTBins, 0.0, dMaxDeltaVertexT);

816

817

// Max Track DOCA

818

locHistName = "MaxTrackDOCA";

819

locHistTitle = locStepROOTName + string(";Largest Track DOCA (cm)");

820

dHistDeque_MaxTrackDOCA[loc_i] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDOCABins, 0.0, dMaxDOCA);

821

822

for(size_t loc_j = 0; loc_j < locDetectedChargedPIDs.size(); ++loc_j)

823

{

824

locPID = locDetectedChargedPIDs[loc_j];

825

if(dHistDeque_TrackZToCommon[loc_i].find(locPID) != dHistDeque_TrackZToCommon[loc_i].end())

826

continue; //already created for this pid

827

locParticleName = ParticleType(locPID);

828

locParticleROOTName = ParticleName_ROOT(locPID);

829

830

// TrackZ To Common

831

locHistName = string("TrackZToCommon_") + locParticleName;

832

locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";#DeltaVertex-Z (Track, Common) (cm)");

833

dHistDeque_TrackZToCommon[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaVertexZBins, dMinDeltaVertexZ, dMaxDeltaVertexZ);

834

835

// TrackT To Common

836

locHistName = string("TrackTToCommon_") + locParticleName;

837

locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";#DeltaVertex-T (Track, Common) (ns)");

838

dHistDeque_TrackTToCommon[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaVertexTBins, dMinDeltaVertexT, dMaxDeltaVertexT);

839

840

// TrackDOCA To Common

841

locHistName = string("TrackDOCAToCommon_") + locParticleName;

842

locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";DOCA (Track, Common) (cm)");

843

dHistDeque_TrackDOCAToCommon[loc_i][locPID] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDOCABins, dMinDOCA, dMaxDOCA);

844

845

// DeltaT Vs P against beam photon

846

if((loc_i == 0) && DAnalysis::Get_IsFirstStepBeam(Get_Reaction()) && (dHistMap_BeamTrackDeltaTVsP.find(locPID) == dHistMap_BeamTrackDeltaTVsP.end()))

847

{

848

locHistName = string("TrackDeltaTVsP_") + ParticleType(locPID) + string("_Beam") + ParticleType(locPID);

849

locHistTitle = locStepROOTName + string(";") + ParticleName_ROOT(locPID) + string(" Momentum (GeV/c);t_{") + ParticleName_ROOT(locPID) + string("} - t_{Beam ") + ParticleName_ROOT(locPID) + string("} (ns)");

850

dHistMap_BeamTrackDeltaTVsP[locPID] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumDeltaVertexTBins, dMinDeltaVertexT, dMaxDeltaVertexT);

851

}

852

}

853

854

//delta-t vs p

855

auto locDetectedChargedPIDs_HasDupes = Get_Reaction()->Get_FinalPIDs(loc_i, false, false, d_Charged, true);

856

for(int loc_j = 0; loc_j < int(locDetectedChargedPIDs_HasDupes.size()) - 1; ++loc_j)

857

{

858

locPID = locDetectedChargedPIDs_HasDupes[loc_j];

859

for(size_t loc_k = loc_j + 1; loc_k < locDetectedChargedPIDs_HasDupes.size(); ++loc_k)

860

{

861

if(ParticleMass(locDetectedChargedPIDs_HasDupes[loc_k]) > ParticleMass(locPID))

862

{

863

locHigherMassPID = locDetectedChargedPIDs_HasDupes[loc_k];

864

locLowerMassPID = locPID;

865

}

866

else

867

{

868

locHigherMassPID = locPID;

869

locLowerMassPID = locDetectedChargedPIDs_HasDupes[loc_k];

870

}

871

pair<Particle_t, Particle_t> locParticlePair(locHigherMassPID, locLowerMassPID);

872

if(dHistDeque_TrackDeltaTVsP[loc_i].find(locParticlePair) != dHistDeque_TrackDeltaTVsP[loc_i].end())

873

continue; //already created for this pair

874

875

locHigherMassParticleName = ParticleType(locHigherMassPID);

876

locHigherMassParticleROOTName = ParticleName_ROOT(locHigherMassPID);

877

locLowerMassParticleName = ParticleType(locLowerMassPID);

878

locLowerMassParticleROOTName = ParticleName_ROOT(locLowerMassPID);

879

880

// DeltaT Vs P

881

locHistName = string("TrackDeltaTVsP_") + locHigherMassParticleName + string("_") + locLowerMassParticleName;

882

locHistTitle = locStepROOTName + string(";") + locHigherMassParticleROOTName + string(" Momentum (GeV/c);t_{") + locHigherMassParticleROOTName + string("} - t_{") + locLowerMassParticleROOTName + string("} (ns)");

883

dHistDeque_TrackDeltaTVsP[loc_i][locParticlePair] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumDeltaVertexTBins, dMinDeltaVertexT, dMaxDeltaVertexT);

884

}

885

}

886

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

887

}

888

//Return to the base directory

889

ChangeTo_BaseDirectory();

890

}

891

japp->RootUnLock(); //RELEASE ROOT LOCK!!

892

}

893

894

void DHistogramAction_TrackVertexComparison::Run_Update(JEventLoop* locEventLoop)

895

{

896

vector<const DAnalysisUtilities*> locAnalysisUtilitiesVector;

897

locEventLoop->Get(locAnalysisUtilitiesVector);

898

dAnalysisUtilities = locAnalysisUtilitiesVector[0];

899

}

900

901

bool DHistogramAction_TrackVertexComparison::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)

902

{

903

Particle_t locPID;

904

double locDOCA, locDeltaZ, locDeltaT, locMaxDOCA, locMaxDeltaZ, locMaxDeltaT;

905

906

//note: very difficult to tell when results will be duplicate: just histogram all combos

907

for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i)

908

{

909

const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(loc_i);

910

auto locParticles = locParticleComboStep->Get_FinalParticles_Measured(Get_Reaction()->Get_ReactionStep(loc_i), d_Charged);

911

if(locParticles.empty())

912

continue;

913

914

//Grab/Find common vertex & time

915

DVector3 locVertex = locParticleComboStep->Get_Position();

916

double locVertexTime = locParticleComboStep->Get_Time();

917

918

for(size_t loc_j = 0; loc_j < locParticles.size(); ++loc_j)

919

{

920

locPID = locParticles[loc_j]->PID();

921

922

//find max's

923

locMaxDOCA = -1.0;

924

locMaxDeltaZ = -1.0;

925

locMaxDeltaT = -1.0;

926

for(size_t loc_k = loc_j + 1; loc_k < locParticles.size(); ++loc_k)

927

{

928

locDeltaZ = fabs(locParticles[loc_j]->position().Z() - locParticles[loc_k]->position().Z());

929

if(locDeltaZ > locMaxDeltaZ)

930

locMaxDeltaZ = locDeltaZ;

931

932

locDeltaT = fabs(locParticles[loc_j]->time() - locParticles[loc_k]->time());

933

if(locDeltaT > locMaxDeltaT)

934

locMaxDeltaT = locDeltaT;

935

936

locDOCA = dAnalysisUtilities->Calc_DOCA(locParticles[loc_j], locParticles[loc_k]);

937

if(locDOCA > locMaxDOCA)

938

locMaxDOCA = locDOCA;

939

}

940

941

//delta-t vs p

942

deque<pair<const DKinematicData*, size_t> > locParticlePairs;

943

size_t locHigherMassParticleIndex, locLowerMassParticleIndex;

944

//minimize number of locks by keeping track of the results and saving them at the end

945

deque<pair<Particle_t, Particle_t> > locPIDPairs;

946

deque<double> locPs;

947

deque<double> locDeltaTs;

948

for(size_t loc_k = loc_j + 1; loc_k < locParticles.size(); ++loc_k)

949

{

950

locParticlePairs.clear();

951

locParticlePairs.push_back(pair<const DKinematicData*, size_t>(locParticles[loc_j], loc_i));

952

locParticlePairs.push_back(pair<const DKinematicData*, size_t>(locParticles[loc_k], loc_i));

953

954

if(locParticles[loc_k]->mass() > locParticles[loc_j]->mass())

955

{

956

locHigherMassParticleIndex = loc_k;

957

locLowerMassParticleIndex = loc_j;

958

}

959

else

960

{

961

locHigherMassParticleIndex = loc_j;

962

locLowerMassParticleIndex = loc_k;

963

}

964

965

locDeltaTs.push_back(locParticles[locHigherMassParticleIndex]->time() - locParticles[locLowerMassParticleIndex]->time());

966

locPs.push_back(locParticles[locHigherMassParticleIndex]->momentum().Mag());

967

locPIDPairs.push_back(pair<Particle_t, Particle_t>(locParticles[locHigherMassParticleIndex]->PID(), locParticles[locLowerMassParticleIndex]->PID()));

968

}

969

970

const DKinematicData* locBeamParticle = ((loc_i == 0) && DAnalysis::Get_IsFirstStepBeam(Get_Reaction())) ? locParticleComboStep->Get_InitialParticle() : NULL__null;

971

double locBeamDeltaT = (locBeamParticle != NULL__null) ? locParticles[loc_j]->time() - locBeamParticle->time() : numeric_limits<double>::quiet_NaN();

972

locDOCA = dAnalysisUtilities->Calc_DOCAToVertex(locParticles[loc_j], locVertex);

973

974

//do all at once to reduce #locks & amount of time within the lock

975

//FILL HISTOGRAMS

976

//Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock

977

//Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex

978

Lock_Action();

979

{

980

//comparison to common vertex/time

981

dHistDeque_TrackZToCommon[loc_i][locPID]->Fill(locParticles[loc_j]->position().Z() - locVertex.Z());

982

dHistDeque_TrackTToCommon[loc_i][locPID]->Fill(locParticles[loc_j]->time() - locVertexTime);

983

dHistDeque_TrackDOCAToCommon[loc_i][locPID]->Fill(locDOCA);

984

//hist max's

985

if(locMaxDeltaZ > 0.0) //else none found (e.g. only 1 detected charged track)

986

{

987

dHistDeque_MaxTrackDeltaZ[loc_i]->Fill(locMaxDeltaZ);

988

dHistDeque_MaxTrackDeltaT[loc_i]->Fill(locMaxDeltaT);

989

dHistDeque_MaxTrackDOCA[loc_i]->Fill(locMaxDOCA);

990

}

991

//delta-t's

992

if(locBeamParticle != NULL__null)

993

dHistMap_BeamTrackDeltaTVsP[locPID]->Fill(locParticles[loc_j]->momentum().Mag(), locBeamDeltaT);

994

for(size_t loc_k = 0; loc_k < locPIDPairs.size(); ++loc_k)

995

{

996

if(dHistDeque_TrackDeltaTVsP[loc_i].find(locPIDPairs[loc_k]) == dHistDeque_TrackDeltaTVsP[loc_i].end())

997

{

998

//pair not found: equal masses and order switched somehow //e.g. mass set differently between REST and reconstruction

999

pair<Particle_t, Particle_t> locTempPIDPair(locPIDPairs[loc_k]);

1000

locPIDPairs[loc_k].first = locTempPIDPair.second;

1001

locPIDPairs[loc_k].second = locTempPIDPair.first;

1002

}

1003

dHistDeque_TrackDeltaTVsP[loc_i][locPIDPairs[loc_k]]->Fill(locPs[loc_k], locDeltaTs[loc_k]);

1004

}

1005

}

1006

Unlock_Action();

1007

} //end of particle loop

1008

} //end of step loop

1009

return true;

1010

}

1011

1012

void DHistogramAction_ParticleComboKinematics::Initialize(JEventLoop* locEventLoop)

1013

{

1014

if(Get_UseKinFitResultsFlag() && (Get_Reaction()->Get_KinFitType() == d_NoFit))

1015

{

1016

cout << "WARNING: REQUESTED HISTOGRAM OF KINEMATIC FIT RESULTS WHEN NO KINEMATIC FIT!!!" << endl;

1017

return; //no fit performed, but kinfit data requested!!

1018

}

1019

1020

string locHistName, locHistTitle, locStepROOTName, locParticleName, locParticleROOTName;

1021

Particle_t locPID;

1022

1023

Run_Update(locEventLoop);

1024

1025

size_t locNumSteps = Get_Reaction()->Get_NumReactionSteps();

1026

dHistDeque_PVsTheta.resize(locNumSteps);

1027

dHistDeque_PhiVsTheta.resize(locNumSteps);

1028

dHistDeque_BetaVsP.resize(locNumSteps);

1029

dHistDeque_DeltaBetaVsP.resize(locNumSteps);

1030

dHistDeque_P.resize(locNumSteps);

1031

dHistDeque_Theta.resize(locNumSteps);

1032

dHistDeque_Phi.resize(locNumSteps);

1033

dHistDeque_VertexZ.resize(locNumSteps);

1034

dHistDeque_VertexYVsX.resize(locNumSteps);

1035

1036

auto locIsFirstStepBeam = DAnalysis::Get_IsFirstStepBeam(Get_Reaction());

1037

1038

//CREATE THE HISTOGRAMS

1039

//Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock

1040

japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!

1041

{

1042

CreateAndChangeTo_ActionDirectory();

1043

1044

//beam particle

1045

locPID = Get_Reaction()->Get_ReactionStep(0)->Get_InitialPID();

1046

if(locIsFirstStepBeam)

1047

{

1048

locParticleName = string("Beam_") + ParticleType(locPID);

1049

CreateAndChangeTo_Directory(locParticleName, locParticleName);

1050

locParticleROOTName = ParticleName_ROOT(locPID);

1051

1052

// Momentum

1053

locHistName = "Momentum";

1054

locHistTitle = string("Beam ") + locParticleROOTName + string(";p (GeV/c)");

1055

dBeamParticleHist_P = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPBins, dMinP, dMaxP);

1056

1057

// Theta

1058

locHistName = "Theta";

1059

locHistTitle = string("Beam ") + locParticleROOTName + string(";#theta#circ");

1060

dBeamParticleHist_Theta = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumThetaBins, dMinTheta, dMaxTheta);

1061

1062

// Phi

1063

locHistName = "Phi";

1064

locHistTitle = string("Beam ") + locParticleROOTName + string(";#phi#circ");

1065

dBeamParticleHist_Phi = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPhiBins, dMinPhi, dMaxPhi);

1066

1067

// P Vs Theta

1068

locHistName = "PVsTheta";

1069

locHistTitle = string("Beam ") + locParticleROOTName + string(";#theta#circ;p (GeV/c)");

1070

dBeamParticleHist_PVsTheta = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);

1071

1072

// Phi Vs Theta

1073

locHistName = "PhiVsTheta";

1074

locHistTitle = string("Beam ") + locParticleROOTName + string(";#theta#circ;#phi#circ");

1075

dBeamParticleHist_PhiVsTheta = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPhiBins, dMinPhi, dMaxPhi);

1076

1077

// Vertex-Z

1078

locHistName = "VertexZ";

1079

locHistTitle = string("Beam ") + locParticleROOTName + string(";Vertex-Z (cm)");

1080

dBeamParticleHist_VertexZ = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumVertexZBins, dMinVertexZ, dMaxVertexZ);

1081

1082

// Vertex-Y Vs Vertex-X

1083

locHistName = "VertexYVsX";

1084

locHistTitle = string("Beam ") + locParticleROOTName + string(";Vertex-X (cm);Vertex-Y (cm)");

1085

dBeamParticleHist_VertexYVsX = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY);

1086

1087

// Delta-T (Beam, RF)

1088

locHistName = "DeltaTRF";

1089

locHistTitle = string("Beam ") + locParticleROOTName + string(";#Deltat_{Beam - RF} (ns)");

1090

dBeamParticleHist_DeltaTRF = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumDeltaTRFBins, dMinDeltaTRF, dMaxDeltaTRF);

1091

1092

// Delta-T (Beam, RF) Vs Beam E

1093

locHistName = "DeltaTRFVsBeamE";

1094

locHistTitle = string("Beam ") + locParticleROOTName + string(";E (GeV);#Deltat_{Beam - RF} (ns)");

1095

dBeamParticleHist_DeltaTRFVsBeamE = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumDeltaTRFBins, dMinDeltaTRF, dMaxDeltaTRF);

1096

1097

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

1098

}

1099

1100

//Steps

1101

for(size_t loc_i = 0; loc_i < locNumSteps; ++loc_i)

1102

{

1103

const DReactionStep* locReactionStep = Get_Reaction()->Get_ReactionStep(loc_i);

1104

ostringstream locStepName;

1105

locStepName << "Step" << loc_i << "__" << DAnalysis::Get_StepName(locReactionStep, true, false);

1106

locStepROOTName = DAnalysis::Get_StepName(locReactionStep, true, true);

1107

1108

Particle_t locInitialPID = locReactionStep->Get_InitialPID();

1109

1110

1111

//get PIDs

1112

vector<Particle_t> locPIDs;

1113

bool locLastPIDMissingFlag = false;

1114

if(!Get_UseKinFitResultsFlag()) //measured, ignore missing & decaying particles (ignore target anyway)

1115

locPIDs = Get_Reaction()->Get_FinalPIDs(loc_i, false, false, d_AllCharges, false);

1116

else //kinematic fit: decaying & missing particles are reconstructed

1117

{

1118

locPIDs = Get_Reaction()->Get_FinalPIDs(loc_i, false, true, d_AllCharges, false);

1119

if((!locIsFirstStepBeam) || (loc_i != 0)) //add decaying parent particle //skip if on beam particle!

1120

locPIDs.insert(locPIDs.begin(), locInitialPID);

1121

1122

Particle_t locMissingPID = locReactionStep->Get_MissingPID();

1123

if(locMissingPID != Unknown)

1124

{

1125

locPIDs.push_back(locMissingPID);

1126

locLastPIDMissingFlag = true;

1127

}

1128

}

1129

1130

bool locDirectoryCreatedFlag = false;

1131

for(size_t loc_j = 0; loc_j < locPIDs.size(); ++loc_j)

1132

{

1133

locPID = locPIDs[loc_j];

1134

bool locIsMissingFlag = false;

1135

if((loc_j == locPIDs.size() - 1) && locLastPIDMissingFlag)

1136

{

1137

locParticleName = string("(") + ParticleType(locPID) + string(")");

1138

locParticleROOTName = string("(") + ParticleName_ROOT(locPID) + string(")");

1139

locIsMissingFlag = true;

1140

}

1141

else

1142

{

1143

locParticleName = ParticleType(locPID);

1144

locParticleROOTName = ParticleName_ROOT(locPID);

1145

if(dHistDeque_P[loc_i].find(locPID) != dHistDeque_P[loc_i].end())

1146

continue; //pid already done

1147

}

1148

1149

if(!locDirectoryCreatedFlag)

1150

{

1151

CreateAndChangeTo_Directory(locStepName.str(), locStepName.str());

1152

locDirectoryCreatedFlag = true;

1153

}

1154

1155

CreateAndChangeTo_Directory(locParticleName, locParticleName);

1156

1157

// Momentum

1158

locHistName = "Momentum";

1159

locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";p (GeV/c)");

1160

dHistDeque_P[loc_i][locPID][locIsMissingFlag] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPBins, dMinP, dMaxP);

1161

1162

// Theta

1163

locHistName = "Theta";

1164

locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";#theta#circ");

1165

dHistDeque_Theta[loc_i][locPID][locIsMissingFlag] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumThetaBins, dMinTheta, dMaxTheta);

1166

1167

// Phi

1168

locHistName = "Phi";

1169

locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";#phi#circ");

1170

dHistDeque_Phi[loc_i][locPID][locIsMissingFlag] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPhiBins, dMinPhi, dMaxPhi);

1171

1172

// P Vs Theta

1173

locHistName = "PVsTheta";

1174

locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";#theta#circ;p (GeV/c)");

1175

dHistDeque_PVsTheta[loc_i][locPID][locIsMissingFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPBins, dMinP, dMaxP);

1176

1177

// Phi Vs Theta

1178

locHistName = "PhiVsTheta";

1179

locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";#theta#circ;#phi#circ");

1180

dHistDeque_PhiVsTheta[loc_i][locPID][locIsMissingFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DPhiBins, dMinPhi, dMaxPhi);

1181

1182

//beta vs p

1183

locHistName = "BetaVsP";

1184

locHistTitle = locParticleROOTName + string(";p (GeV/c);#beta");

1185

dHistDeque_BetaVsP[loc_i][locPID][locIsMissingFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumBetaBins, dMinBeta, dMaxBeta);

1186

1187

//delta-beta vs p

1188

locHistName = "DeltaBetaVsP";

1189

locHistTitle = locParticleROOTName + string(";p (GeV/c);#Delta#beta");

1190

dHistDeque_DeltaBetaVsP[loc_i][locPID][locIsMissingFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPBins, dMinP, dMaxP, dNumDeltaBetaBins, dMinDeltaBeta, dMaxDeltaBeta);

1191

1192

// Vertex-Z

1193

locHistName = "VertexZ";

1194

locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";Vertex-Z (cm)");

1195

dHistDeque_VertexZ[loc_i][locPID][locIsMissingFlag] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumVertexZBins, dMinVertexZ, dMaxVertexZ);

1196

1197

// Vertex-Y Vs Vertex-X

1198

locHistName = "VertexYVsX";

1199

locHistTitle = locParticleROOTName + string(", ") + locStepROOTName + string(";Vertex-X (cm);Vertex-Y (cm)");

1200

dHistDeque_VertexYVsX[loc_i][locPID][locIsMissingFlag] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY);

1201

1202

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

1203

} //end of particle loop

1204

1205

// Vertex

1206

string locInitParticleROOTName = ParticleName_ROOT(locInitialPID);

1207

string locInitParticleName = ParticleType(locInitialPID);

1208

if((loc_i == 0) || IsDetachedVertex(locInitialPID))

1209

{

1210

if(!locDirectoryCreatedFlag)

1211

{

1212

CreateAndChangeTo_Directory(locStepName.str(), locStepName.str());

1213

locDirectoryCreatedFlag = true;

1214

}

1215

1216

locHistName = "StepVertexZ";

1217

locHistTitle = (loc_i == 0) ? ";Production Vertex-Z (cm)" : string(";") + locInitParticleROOTName + string(" Decay Vertex-Z (cm)");

1218

dHistMap_StepVertexZ[loc_i] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumVertexZBins, dMinVertexZ, dMaxVertexZ);

1219

1220

locHistName = "StepVertexYVsX";

1221

locHistTitle = (loc_i == 0) ? "Production Vertex" : locInitParticleROOTName + string(" Decay Vertex");

1222

locHistTitle += string(";Vertex-X (cm);Vertex-Y (cm)");

1223

dHistMap_StepVertexYVsX[loc_i] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY, dNumVertexXYBins, dMinVertexXY, dMaxVertexXY);

1224

}

1225

1226

if((loc_i != 0) && IsDetachedVertex(locInitialPID))

1227

{

1228

if(!locDirectoryCreatedFlag)

1229

{

1230

CreateAndChangeTo_Directory(locStepName.str(), locStepName.str());

1231

locDirectoryCreatedFlag = true;

1232

}

1233

1234

locHistName = locInitParticleName + string("PathLength");

1235

locHistTitle = string(";") + locInitParticleROOTName + string(" Path Length (cm)");

1236

dHistMap_DetachedPathLength[loc_i] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPathLengthBins, 0.0, dMaxPathLength);

1237

1238

locHistName = locInitParticleName + string("Lifetime");

1239

locHistTitle = string(";") + locInitParticleROOTName + string(" Lifetime (ns)");

1240

dHistMap_DetachedLifetime[loc_i] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumLifetimeBins, 0.0, dMaxLifetime);

1241

1242

if(Get_UseKinFitResultsFlag())

1243

{

1244

locHistName = locInitParticleName + string("Lifetime_RestFrame");

1245

locHistTitle = string(";") + locInitParticleROOTName + string(" Rest Frame Lifetime (ns)");

1246

dHistMap_DetachedLifetime_RestFrame[loc_i] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumLifetimeBins, 0.0, dMaxLifetime);

1247

}

1248

}

1249

1250

if(locDirectoryCreatedFlag)

1251

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

1252

} //end of step loop

1253

1254

//Return to the base directory

1255

ChangeTo_BaseDirectory();

1256

}

1257

japp->RootUnLock(); //RELEASE ROOT LOCK!!

1258

}

1259

1260

void DHistogramAction_ParticleComboKinematics::Run_Update(JEventLoop* locEventLoop)

1261

{

1262

vector<const DParticleID*> locParticleIDs;

1263

locEventLoop->Get(locParticleIDs);

1264

1265

DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());

1266

vector<const DAnalysisUtilities*> locAnalysisUtilitiesVector;

1267

locEventLoop->Get(locAnalysisUtilitiesVector);

1268

1269

DGeometry *locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());

1270

locGeometry->GetTargetZ(dTargetZCenter);

1271

1272

dParticleID = locParticleIDs[0];

1273

dAnalysisUtilities = locAnalysisUtilitiesVector[0];

1274

}

1275

1276

bool DHistogramAction_ParticleComboKinematics::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)

1277

{

1278

if(Get_UseKinFitResultsFlag() && (Get_Reaction()->Get_KinFitType() == d_NoFit))

1279

{

1280

cout << "WARNING: REQUESTED HISTOGRAM OF KINEMATIC FIT RESULTS WHEN NO KINEMATIC FIT!!! Skipping histogram." << endl;

1281

return true; //no fit performed, but kinfit data requested!!

1282

}

1283

1284

const DEventRFBunch* locEventRFBunch = locParticleCombo->Get_EventRFBunch();

1285

auto locIsFirstStepBeam = DAnalysis::Get_IsFirstStepBeam(Get_Reaction());

1286

1287

const DKinematicData* locKinematicData;

1288

for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i)

1289

{

1290

const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(loc_i);

1291

auto locReactionStep = Get_Reaction()->Get_ReactionStep(loc_i);

1292

1293

//initial particle

1294

Particle_t locInitialPID = locReactionStep->Get_InitialPID();

1295

if(Get_UseKinFitResultsFlag())

1296

locKinematicData = locParticleComboStep->Get_InitialParticle();

1297

else

1298

locKinematicData = locParticleComboStep->Get_InitialParticle_Measured();

1299

if(locKinematicData != NULL__null)

1300

{

1301

if(locIsFirstStepBeam && (loc_i == 0))

1302

{

1303

//check if will be duplicate

1304

const JObject* locSourceObject = locParticleComboStep->Get_InitialParticle_Measured();

1305

if(dPreviouslyHistogrammedBeamParticles.find(locSourceObject) == dPreviouslyHistogrammedBeamParticles.end())

1306

{

1307

dPreviouslyHistogrammedBeamParticles.insert(locSourceObject);

1308

Fill_BeamHists(locKinematicData, locEventRFBunch);

1309

}

1310

}

1311

else if(Get_UseKinFitResultsFlag()) //decaying particle, but kinfit so can hist

1312

Fill_Hists(locEventLoop, locKinematicData, false, loc_i); //has many source object, and is unique to this combo: no dupes to check against: let it ride

1313

}

1314

1315

//VERTEX INFORMATION

1316

//other than first, skipped if not detached vertex

1317

DLorentzVector locStepSpacetimeVertex = locParticleComboStep->Get_SpacetimeVertex();

1318

if((loc_i == 0) || IsDetachedVertex(locInitialPID))

1319

{

1320

dHistMap_StepVertexZ[loc_i]->Fill(locStepSpacetimeVertex.Z());

1321

dHistMap_StepVertexYVsX[loc_i]->Fill(locStepSpacetimeVertex.X(), locStepSpacetimeVertex.Y());

1322

}

1323

1324

//DETACHED VERTEX INFORMATION

1325

if((loc_i != 0) && IsDetachedVertex(locInitialPID))

1326

{

1327

int locFromStepIndex = DAnalysis::Get_InitialParticleDecayFromIndices(Get_Reaction(), loc_i).first;

1328

DLorentzVector locFromSpacetimeVertex = locParticleCombo->Get_ParticleComboStep(locFromStepIndex)->Get_SpacetimeVertex();

1329

DLorentzVector locDeltaSpacetime = locStepSpacetimeVertex - locFromSpacetimeVertex;

1330

1331

double locPathLength = locDeltaSpacetime.Vect().Mag();

1332

dHistMap_DetachedPathLength[loc_i]->Fill(locPathLength);

1333

dHistMap_DetachedLifetime[loc_i]->Fill(locDeltaSpacetime.T());

1334

1335

if(locKinematicData != NULL__null)

1336

{

1337

DLorentzVector locInitialP4 = locKinematicData->lorentzMomentum();

1338

//below, t, x, and tau are really delta-t, delta-x, and delta-tau

1339

//tau = gamma*(t - beta*x/c) //plug in: t = x/(beta*c), gamma = 1/sqrt(1 - beta^2)

1340

//tau = (x/(beta*c) - beta*x/c)/sqrt(1 - beta^2) //plug in beta = p*c/E

1341

//tau = (x*E/(p*c^2) - p*x/E)/sqrt(1 - p^2*c^2/E^2) //multiply num & denom by E*P, factor x/c^2 out of num

1342

//tau = (x/c^2)*(E^2 - p^2*c^2)/(p*sqrt(E^2 - p^2*c^2)) //plug in m^2*c^4 = E^2 - p^2*c^2

1343

//tau = (x/c^2)*m^2*c^4/(p*m*c^2) //cancel c's & m's

1344

//tau = x*m/p

1345

//however, in data, p & m are in units with c = 1, so need an extra 1/c: tau = x*m/(c*p)

1346

double locRestFrameLifetime = locPathLength*ParticleMass(locInitialPID)/(29.9792458*locInitialP4.P()); //tau

1347

dHistMap_DetachedLifetime_RestFrame[loc_i]->Fill(locRestFrameLifetime);

1348

//note that tau = hbar / Gamma, hbar = 6.582119E-22 MeV*s, Gamma = Resonance FWHM

1349

}

1350

}

1351

1352

//final particles

1353

for(size_t loc_j = 0; loc_j < locParticleComboStep->Get_NumFinalParticles(); ++loc_j)

1354

{

1355

if(Get_UseKinFitResultsFlag())

1356

locKinematicData = locParticleComboStep->Get_FinalParticle(loc_j);

1357

else

1358

locKinematicData = locParticleComboStep->Get_FinalParticle_Measured(loc_j);

1359

if(locKinematicData == NULL__null)

1360

continue; //e.g. a decaying or missing particle whose params aren't set yet

1361

1362

//check if duplicate

1363

const JObject* locSourceObject = Get_FinalParticle_SourceObject(locKinematicData);

1364

if(locSourceObject != NULL__null) //else is reconstructed missing/decaying particle: has many source object, and is unique to this combo: no dupes to check against: let it ride

1365

{

1366

pair<Particle_t, const JObject*> locParticleInfo(locKinematicData->PID(), locSourceObject);

1367

pair<size_t, pair<Particle_t, const JObject*> > locHistInfo(loc_i, locParticleInfo);

1368

if(dPreviouslyHistogrammedParticles.find(locHistInfo) != dPreviouslyHistogrammedParticles.end())

1369

continue; //previously histogrammed

1370

dPreviouslyHistogrammedParticles.insert(locHistInfo);

1371

}

1372

1373

bool locIsMissingFlag = (Get_Reaction()->Get_ReactionStep(loc_i)->Get_MissingParticleIndex() == int(loc_j));

1374

Fill_Hists(locEventLoop, locKinematicData, locIsMissingFlag, loc_i);

1375

} //end of particle loop

1376

} //end of step loop

1377

return true;

1378

}

1379

1380

void DHistogramAction_ParticleComboKinematics::Fill_Hists(JEventLoop* locEventLoop, const DKinematicData* locKinematicData, bool locIsMissingFlag, size_t locStepIndex)

1381

{

1382

Particle_t locPID = locKinematicData->PID();

1383

DVector3 locMomentum = locKinematicData->momentum();

1384

DVector3 locPosition = locKinematicData->position();

1385

1386

double locPhi = locMomentum.Phi()*180.0/TMath::Pi();

1387

double locTheta = locMomentum.Theta()*180.0/TMath::Pi();

1388

double locP = locMomentum.Mag();

1389

1390

double locBeta_Timing = 0.0;

1391

auto locNeutralHypo = dynamic_cast<const DNeutralParticleHypothesis*>(locKinematicData);

1392

if(locNeutralHypo != nullptr)

1393

locBeta_Timing = locNeutralHypo->measuredBeta();

1394

else

1395

{

1396

auto locChargedHypo = dynamic_cast<const DChargedTrackHypothesis*>(locKinematicData);

1397

if(locChargedHypo != nullptr)

1398

locBeta_Timing = locChargedHypo->measuredBeta();

1399

}

1400

double locDeltaBeta = locBeta_Timing - locKinematicData->lorentzMomentum().Beta();

1401

1402

//FILL HISTOGRAMS

1403

//Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock

1404

//Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex

1405

Lock_Action();

1406

{

1407

dHistDeque_P[locStepIndex][locPID][locIsMissingFlag]->Fill(locP);

1408

dHistDeque_Phi[locStepIndex][locPID][locIsMissingFlag]->Fill(locPhi);

1409

dHistDeque_Theta[locStepIndex][locPID][locIsMissingFlag]->Fill(locTheta);

1410

dHistDeque_PVsTheta[locStepIndex][locPID][locIsMissingFlag]->Fill(locTheta, locP);

1411

dHistDeque_PhiVsTheta[locStepIndex][locPID][locIsMissingFlag]->Fill(locTheta, locPhi);

1412

dHistDeque_BetaVsP[locStepIndex][locPID][locIsMissingFlag]->Fill(locP, locBeta_Timing);

1413

dHistDeque_DeltaBetaVsP[locStepIndex][locPID][locIsMissingFlag]->Fill(locP, locDeltaBeta);

1414

dHistDeque_VertexZ[locStepIndex][locPID][locIsMissingFlag]->Fill(locKinematicData->position().Z());

1415

dHistDeque_VertexYVsX[locStepIndex][locPID][locIsMissingFlag]->Fill(locKinematicData->position().X(), locKinematicData->position().Y());

1416

}

1417

Unlock_Action();

1418

}

1419

1420

void DHistogramAction_ParticleComboKinematics::Fill_BeamHists(const DKinematicData* locKinematicData, const DEventRFBunch* locEventRFBunch)

1421

{

1422

DVector3 locMomentum = locKinematicData->momentum();

1423

double locPhi = locMomentum.Phi()*180.0/TMath::Pi();

1424

double locTheta = locMomentum.Theta()*180.0/TMath::Pi();

1425

double locP = locMomentum.Mag();

1426

double locDeltaTRF = locKinematicData->time() - (locEventRFBunch->dTime + (locKinematicData->z() - dTargetZCenter)/29.9792458);

1427

1428

//FILL HISTOGRAMS

1429

//Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock

1430

//Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex

1431

Lock_Action();

1432

{

1433

dBeamParticleHist_P->Fill(locP);

1434

dBeamParticleHist_Phi->Fill(locPhi);

1435

dBeamParticleHist_Theta->Fill(locTheta);

1436

dBeamParticleHist_PVsTheta->Fill(locTheta, locP);

1437

dBeamParticleHist_PhiVsTheta->Fill(locTheta, locPhi);

1438

dBeamParticleHist_VertexZ->Fill(locKinematicData->position().Z());

1439

dBeamParticleHist_VertexYVsX->Fill(locKinematicData->position().X(), locKinematicData->position().Y());

1440

dBeamParticleHist_DeltaTRF->Fill(locDeltaTRF);

1441

dBeamParticleHist_DeltaTRFVsBeamE->Fill(locKinematicData->energy(), locDeltaTRF);

1442

}

1443

Unlock_Action();

1444

}

1445

1446

void DHistogramAction_InvariantMass::Initialize(JEventLoop* locEventLoop)

1447

{

1448

string locHistName, locHistTitle;

1449

double locMassPerBin = 1000.0*(dMaxMass - dMinMass)/((double)dNumMassBins);

1450

1451

Run_Update(locEventLoop);

1452

1453

string locParticleNamesForHist = "";

1454

if(dInitialPID != Unknown)

1455

{

1456

auto locChainPIDs = DAnalysis::Get_ChainPIDs(Get_Reaction(), dInitialPID, !Get_UseKinFitResultsFlag(), true);

1457

locParticleNamesForHist = DAnalysis::Convert_PIDsToROOTName(locChainPIDs);

1458

}

1459

else

1460

{

1461

for(size_t loc_i = 0; loc_i < dToIncludePIDs.size(); ++loc_i)

1462

locParticleNamesForHist += ParticleName_ROOT(dToIncludePIDs[loc_i]);

1463

}

1464

1465

bool locBeamPresent = DAnalysis::Get_IsFirstStepBeam(Get_Reaction());

1466

1467

//CREATE THE HISTOGRAMS

1468

//Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock

1469

japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!

1470

{

1471

CreateAndChangeTo_ActionDirectory();

1472

1473

locHistName = "InvariantMass";

1474

ostringstream locStream;

1475

locStream << locMassPerBin;

1476

locHistTitle = string(";") + locParticleNamesForHist + string(" Invariant Mass (GeV/c^{2});# Combos / ") + locStream.str() + string(" MeV/c^{2}");

1477

dHist_InvariantMass = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumMassBins, dMinMass, dMaxMass);

1478

1479

if(locBeamPresent)

1480

{

1481

locHistName = "InvariantMassVsBeamE";

1482

locHistTitle = string(";Beam Energy (GeV);") + locParticleNamesForHist + string(" Invariant Mass (GeV/c^{2})");

1483

dHist_InvariantMassVsBeamE = GetOrCreate_Histogram<TH2D>(locHistName, locHistTitle, dNum2DBeamEBins, dMinBeamE, dMaxBeamE, dNum2DMassBins, dMinMass, dMaxMass);

1484

}

1485

1486

//Return to the base directory

1487

ChangeTo_BaseDirectory();

1488

}

1489

japp->RootUnLock(); //RELEASE ROOT LOCK!!

1490

}

1491

1492

void DHistogramAction_InvariantMass::Run_Update(JEventLoop* locEventLoop)

1493

{

1494

DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());

1495

DGeometry *locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());

1496

locGeometry->GetTargetZ(dTargetZCenter);

1497

1498

locEventLoop->GetSingle(dAnalysisUtilities);

1499

1500

//BEAM BUNCH PERIOD

1501

vector<double> locBeamPeriodVector;

1502

locEventLoop->GetCalib("PHOTON_BEAM/RF/beam_period", locBeamPeriodVector);

1503

dBeamBunchPeriod = locBeamPeriodVector[0];

1504

}

1505

1506

bool DHistogramAction_InvariantMass::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)

1507

{

1508

auto locFirstStep = locParticleCombo->Get_ParticleComboStep(0);

1509

bool locBeamPresent = DAnalysis::Get_IsFirstStepBeam(Get_Reaction());

1510

auto locBeam = Get_UseKinFitResultsFlag() ? locFirstStep->Get_InitialParticle() : locFirstStep->Get_InitialParticle_Measured();

1511

1512

vector<double> locMassesToFill;

1513

vector<double> loc2DMassesToFill;

1514

for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i)

1515

{

1516

const DReactionStep* locReactionStep = Get_Reaction()->Get_ReactionStep(loc_i);

1517

if((dInitialPID != Unknown) && (locReactionStep->Get_InitialPID() != dInitialPID))

1518

continue;

1519

if((dStepIndex != -1) && (int(loc_i) != dStepIndex))

1520

continue;

1521

1522

//build all possible combinations of the included pids

1523

set<set<size_t> > locIndexCombos = dAnalysisUtilities->Build_IndexCombos(locReactionStep, dToIncludePIDs);

1524

1525

//loop over them

1526

set<set<size_t> >::iterator locComboIterator = locIndexCombos.begin();

1527

for(; locComboIterator != locIndexCombos.end(); ++locComboIterator)

1528

{

1529

set<pair<const JObject*, unsigned int> > locSourceObjects;

1530

DLorentzVector locFinalStateP4 = dAnalysisUtilities->Calc_FinalStateP4(Get_Reaction(), locParticleCombo, loc_i, *locComboIterator, locSourceObjects, Get_UseKinFitResultsFlag());

1531

if(dPreviousSourceObjects.find(locSourceObjects) == dPreviousSourceObjects.end())

1532

{

1533

dPreviousSourceObjects.insert(locSourceObjects);

1534

locMassesToFill.push_back(locFinalStateP4.M());

1535

}

1536

1537

auto locBeamSourceObjects = std::make_pair(locSourceObjects, locBeam);

1538

if(locBeamPresent && (dPreviousSourceObjects_Beam.find(locBeamSourceObjects) == dPreviousSourceObjects_Beam.end()))

1539

{

1540

dPreviousSourceObjects_Beam.insert(locBeamSourceObjects);

1541

loc2DMassesToFill.push_back(locFinalStateP4.M());

1542

}

1543

}

1544

//don't break: e.g. if multiple pi0's, histogram invariant mass of each one

1545

}

1546

1547

// calculate accidental scaling factor (if needed)

1548

double locWeight = 1.;

1549

if(dSubtractAccidentals) {

1550

auto locIsFirstStepBeam = DAnalysis::Get_IsFirstStepBeam(Get_Reaction());

1551

if(locIsFirstStepBeam) {

1552

const DEventRFBunch* locEventRFBunch = locParticleCombo->Get_EventRFBunch();

1553

const DKinematicData* locBeamParticle = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle();

1554

1555

double locDeltaTRF = locBeamParticle->time() - (locEventRFBunch->dTime + (locBeamParticle->z() - dTargetZCenter)/29.9792458);

1556

1557

if(fabs(locDeltaTRF) > dBeamBunchPeriod/2.)

1558

locWeight = -1./(2.*Get_Reaction()->Get_NumPlusMinusRFBunches());

1559

}

1560

}

1561

1562

//FILL HISTOGRAMS

1563

//Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock

1564

//Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex

1565

Lock_Action();

1566

{

1567

if(dSubtractAccidentals) {

1568

for(size_t loc_i = 0; loc_i < locMassesToFill.size(); ++loc_i)

1569

dHist_InvariantMass->Fill(locMassesToFill[loc_i], locWeight);

1570

for(size_t loc_i = 0; loc_i < loc2DMassesToFill.size(); ++loc_i)

1571

dHist_InvariantMassVsBeamE->Fill(locBeam->energy(), loc2DMassesToFill[loc_i], locWeight);

1572

} else {

1573

for(size_t loc_i = 0; loc_i < locMassesToFill.size(); ++loc_i)

1574

dHist_InvariantMass->Fill(locMassesToFill[loc_i]);

1575

for(size_t loc_i = 0; loc_i < loc2DMassesToFill.size(); ++loc_i)

1576

dHist_InvariantMassVsBeamE->Fill(locBeam->energy(), loc2DMassesToFill[loc_i]);

1577

}

1578

}

1579

Unlock_Action();

1580

1581

return true;

1582

}

1583

1584

void DHistogramAction_MissingMass::Initialize(JEventLoop* locEventLoop)

1585

{

1586

double locMassPerBin = 1000.0*(dMaxMass - dMinMass)/((double)dNumMassBins);

1587

string locInitialParticlesROOTName = DAnalysis::Get_InitialParticlesName(Get_Reaction()->Get_ReactionStep(0), true);

1588

vector<Particle_t> locMissingMassOffOfPIDs(dMissingMassOffOfPIDs.begin(), dMissingMassOffOfPIDs.end());

1589

auto locChainPIDs = DAnalysis::Get_ChainPIDs(Get_Reaction(), Get_Reaction()->Get_ReactionStep(0)->Get_InitialPID(), dMissingMassOffOfStepIndex, locMissingMassOffOfPIDs, !Get_UseKinFitResultsFlag(), true);

1590

string locFinalParticlesROOTName = DAnalysis::Convert_PIDsToROOTName(locChainPIDs);

1591

1592

Run_Update(locEventLoop);

1593

1594

//CREATE THE HISTOGRAMS

1595

//Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock

1596

japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!

1597

{

1598

CreateAndChangeTo_ActionDirectory();

1599

1600

string locHistName = "MissingMass";

1601

ostringstream locStream;

1602

locStream << locMassPerBin;

1603

string locHistTitle = string(";") + locInitialParticlesROOTName + string("#rightarrow") + locFinalParticlesROOTName + string(" Missing Mass (GeV/c^{2});# Combos / ") + locStream.str() + string(" MeV/c^{2}");

1604

dHist_MissingMass = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumMassBins, dMinMass, dMaxMass);

1605

1606

locHistName = "MissingMassVsBeamE";

1607

locMassPerBin *= ((double)dNumMassBins)/((double)dNum2DMassBins);

1608

locStream.str("");

1609

locStream << locMassPerBin;

1610

locHistTitle = string(";Beam Energy (GeV);") + locInitialParticlesROOTName + string("#rightarrow") + locFinalParticlesROOTName + string(" Missing Mass (GeV/c^{2})");

1611

dHist_MissingMassVsBeamE = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DBeamEBins, dMinBeamE, dMaxBeamE, dNum2DMassBins, dMinMass, dMaxMass);

1612

1613

locHistName = "MissingMassVsMissingP";

1614

locStream.str("");

1615

locStream << locMassPerBin;

1616

locHistTitle = string(";Missing P (GeV/c);") + locInitialParticlesROOTName + string("#rightarrow") + locFinalParticlesROOTName + string(" Missing Mass (GeV/c^{2})");

1617

dHist_MissingMassVsMissingP = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DMissPBins, dMinMissP, dMaxMissP, dNum2DMassBins, dMinMass, dMaxMass);

1618

1619

//Return to the base directory

1620

ChangeTo_BaseDirectory();

1621

}

1622

japp->RootUnLock(); //RELEASE ROOT LOCK!!

1623

}

1624

1625

void DHistogramAction_MissingMass::Run_Update(JEventLoop* locEventLoop)

1626

{

1627

DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());

1628

DGeometry *locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());

1629

locGeometry->GetTargetZ(dTargetZCenter);

1630

1631

locEventLoop->GetSingle(dAnalysisUtilities);

1632

1633

//BEAM BUNCH PERIOD

1634

vector<double> locBeamPeriodVector;

1635

locEventLoop->GetCalib("PHOTON_BEAM/RF/beam_period", locBeamPeriodVector);

1636

dBeamBunchPeriod = locBeamPeriodVector[0];

1637

}

1638

1639

bool DHistogramAction_MissingMass::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)

1640

{

1641

double locBeamEnergy = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle()->energy();

1642

1643

//build all possible combinations of the included pids

1644

set<set<size_t> > locIndexCombos = dAnalysisUtilities->Build_IndexCombos(Get_Reaction()->Get_ReactionStep(dMissingMassOffOfStepIndex), dMissingMassOffOfPIDs);

1645

1646

//loop over them

1647

vector<pair<double, double> > locMassesToFill; //first is missing mass, second is missing p

1648

set<set<size_t> >::iterator locComboIterator = locIndexCombos.begin();

1649

for(; locComboIterator != locIndexCombos.end(); ++locComboIterator)

1650

{

1651

set<pair<const JObject*, unsigned int> > locSourceObjects;

1652

DLorentzVector locMissingP4 = dAnalysisUtilities->Calc_MissingP4(Get_Reaction(), locParticleCombo, 0, dMissingMassOffOfStepIndex, *locComboIterator, locSourceObjects, Get_UseKinFitResultsFlag());

1653

1654

if(dPreviousSourceObjects.find(locSourceObjects) != dPreviousSourceObjects.end())

1655

continue; //dupe: already histed!

1656

dPreviousSourceObjects.insert(locSourceObjects);

1657

1658

locMassesToFill.push_back(pair<double, double>(locMissingP4.M(), locMissingP4.P()));

1659

}

1660

1661

// calculate accidental scaling factor (if needed)

1662

double locWeight = 1.;

1663

if(dSubtractAccidentals) {

1664

auto locIsFirstStepBeam = DAnalysis::Get_IsFirstStepBeam(Get_Reaction());

1665

if(locIsFirstStepBeam) {

1666

const DEventRFBunch* locEventRFBunch = locParticleCombo->Get_EventRFBunch();

1667

const DKinematicData* locBeamParticle = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle();

1668

1669

double locDeltaTRF = locBeamParticle->time() - (locEventRFBunch->dTime + (locBeamParticle->z() - dTargetZCenter)/29.9792458);

1670

1671

if(fabs(locDeltaTRF) > dBeamBunchPeriod/2.)

1672

locWeight = -1./(2.*Get_Reaction()->Get_NumPlusMinusRFBunches());

1673

}

1674

}

1675

1676

//FILL HISTOGRAMS

1677

//Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock

1678

//Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex

1679

Lock_Action();

1680

{

1681

for(size_t loc_i = 0; loc_i < locMassesToFill.size(); ++loc_i)

1682

{

1683

// only fill histograms with weights if needed, to avoid complications

1684

if(dSubtractAccidentals) {

1685

dHist_MissingMass->Fill(locMassesToFill[loc_i].first, locWeight);

1686

dHist_MissingMassVsBeamE->Fill(locBeamEnergy, locMassesToFill[loc_i].first, locWeight);

1687

dHist_MissingMassVsMissingP->Fill(locMassesToFill[loc_i].second, locMassesToFill[loc_i].first, locWeight);

1688

} else {

1689

dHist_MissingMass->Fill(locMassesToFill[loc_i].first);

1690

dHist_MissingMassVsBeamE->Fill(locBeamEnergy, locMassesToFill[loc_i].first);

1691

dHist_MissingMassVsMissingP->Fill(locMassesToFill[loc_i].second, locMassesToFill[loc_i].first);

1692

}

1693

}

1694

}

1695

Unlock_Action();

1696

1697

return true;

1698

}

1699

1700

void DHistogramAction_MissingMassSquared::Initialize(JEventLoop* locEventLoop)

1701

{

1702

double locMassSqPerBin = 1000.0*1000.0*(dMaxMassSq - dMinMassSq)/((double)dNumMassBins);

1703

string locInitialParticlesROOTName = DAnalysis::Get_InitialParticlesName(Get_Reaction()->Get_ReactionStep(0), true);

1704

vector<Particle_t> locMissingMassOffOfPIDs(dMissingMassOffOfPIDs.begin(), dMissingMassOffOfPIDs.end());

1705

1706

string locFinalParticlesROOTName = DAnalysis::Convert_PIDsToROOTName(locChainPIDs);

1707

1708

Run_Update(locEventLoop);

1709

1710

//CREATE THE HISTOGRAMS

1711

//Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock

1712

japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!

1713

{

1714

CreateAndChangeTo_ActionDirectory();

1715

1716

string locHistName = "MissingMassSquared";

1717

ostringstream locStream;

1718

locStream << locMassSqPerBin;

1719

string locHistTitle = string(";") + locInitialParticlesROOTName + string("#rightarrow") + locFinalParticlesROOTName + string(" Missing Mass Squared (GeV/c^{2})^{2};# Combos / ") + locStream.str() + string(" (MeV/c^{2})^{2}");

1720

dHist_MissingMassSquared = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumMassBins, dMinMassSq, dMaxMassSq);

1721

1722

locHistName = "MissingMassSquaredVsBeamE";

1723

locMassSqPerBin *= ((double)dNumMassBins)/((double)dNum2DMassBins);

1724

locStream.str();

1725

locStream << locMassSqPerBin;

1726

locHistTitle = string(";Beam Energy (GeV);") + locInitialParticlesROOTName + string("#rightarrow") + locFinalParticlesROOTName + string(" Missing Mass Squared (GeV/c^{2})^{2};");

1727

dHist_MissingMassSquaredVsBeamE = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DBeamEBins, dMinBeamE, dMaxBeamE, dNum2DMassBins, dMinMassSq, dMaxMassSq);

1728

1729

locHistName = "MissingMassSquaredVsMissingP";

1730

locStream.str("");

1731

locStream << locMassSqPerBin;

1732

locHistTitle = string(";Missing P (GeV/c);") + locInitialParticlesROOTName + string("#rightarrow") + locFinalParticlesROOTName + string(" Missing Mass Squared (GeV/c^{2})^{2}");

1733

dHist_MissingMassSquaredVsMissingP = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DMissPBins, dMinMissP, dMaxMissP, dNum2DMassBins, dMinMassSq, dMaxMassSq);

1734

1735

//Return to the base directory

1736

ChangeTo_BaseDirectory();

1737

}

1738

japp->RootUnLock(); //RELEASE ROOT LOCK!!

1739

}

1740

1741

void DHistogramAction_MissingMassSquared::Run_Update(JEventLoop* locEventLoop)

1742

{

1743

DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());

1744

DGeometry *locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());

1745

locGeometry->GetTargetZ(dTargetZCenter);

1746

1747

locEventLoop->GetSingle(dAnalysisUtilities);

1748

1749

//BEAM BUNCH PERIOD

1750

vector<double> locBeamPeriodVector;

1751

locEventLoop->GetCalib("PHOTON_BEAM/RF/beam_period", locBeamPeriodVector);

1752

dBeamBunchPeriod = locBeamPeriodVector[0];

1753

}

1754

1755

bool DHistogramAction_MissingMassSquared::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)

1756

{

1757

double locBeamEnergy = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle()->energy();

1758

1759

//build all possible combinations of the included pids

1760

set<set<size_t> > locIndexCombos = dAnalysisUtilities->Build_IndexCombos(Get_Reaction()->Get_ReactionStep(dMissingMassOffOfStepIndex), dMissingMassOffOfPIDs);

1761

1762

//loop over them

1763

vector<pair<double, double> > locMassesToFill; //first is missing mass, second is missing p

1764

set<set<size_t> >::iterator locComboIterator = locIndexCombos.begin();

1765

for(; locComboIterator != locIndexCombos.end(); ++locComboIterator)

1766

{

1767

set<pair<const JObject*, unsigned int> > locSourceObjects;

1768

DLorentzVector locMissingP4 = dAnalysisUtilities->Calc_MissingP4(Get_Reaction(), locParticleCombo, 0, dMissingMassOffOfStepIndex, *locComboIterator, locSourceObjects, Get_UseKinFitResultsFlag());

1769

1770

if(dPreviousSourceObjects.find(locSourceObjects) != dPreviousSourceObjects.end())

1771

continue; //dupe: already histed!

1772

dPreviousSourceObjects.insert(locSourceObjects);

1773

1774

locMassesToFill.push_back(pair<double, double>(locMissingP4.M2(), locMissingP4.P()));

1775

}

1776

1777

// calculate accidental scaling factor (if needed)

1778

double locWeight = 1.;

1779

if(dSubtractAccidentals) {

1780

auto locIsFirstStepBeam = DAnalysis::Get_IsFirstStepBeam(Get_Reaction());

1781

if(locIsFirstStepBeam) {

1782

const DEventRFBunch* locEventRFBunch = locParticleCombo->Get_EventRFBunch();

1783

const DKinematicData* locBeamParticle = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle();

1784

1785

double locDeltaTRF = locBeamParticle->time() - (locEventRFBunch->dTime + (locBeamParticle->z() - dTargetZCenter)/29.9792458);

1786

1787

if(fabs(locDeltaTRF) > dBeamBunchPeriod/2.)

1788

locWeight = -1./(2.*Get_Reaction()->Get_NumPlusMinusRFBunches());

1789

}

1790

}

1791

1792

//FILL HISTOGRAMS

1793

//Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock

1794

//Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex

1795

Lock_Action();

1796

{

1797

for(size_t loc_i = 0; loc_i < locMassesToFill.size(); ++loc_i)

1798

{

1799

if(dSubtractAccidentals) {

1800

dHist_MissingMassSquared->Fill(locMassesToFill[loc_i].first, locWeight);

1801

dHist_MissingMassSquaredVsBeamE->Fill(locBeamEnergy, locMassesToFill[loc_i].first, locWeight);

1802

dHist_MissingMassSquaredVsMissingP->Fill(locMassesToFill[loc_i].second, locMassesToFill[loc_i].first, locWeight);

1803

} else {

1804

dHist_MissingMassSquared->Fill(locMassesToFill[loc_i].first);

1805

dHist_MissingMassSquaredVsBeamE->Fill(locBeamEnergy, locMassesToFill[loc_i].first);

1806

dHist_MissingMassSquaredVsMissingP->Fill(locMassesToFill[loc_i].second, locMassesToFill[loc_i].first);

1807

}

1808

}

1809

}

1810

Unlock_Action();

1811

1812

return true;

1813

}

1814

1815

void DHistogramAction_2DInvariantMass::Initialize(JEventLoop* locEventLoop)

1816

{

1817

string locHistName, locHistTitle;

1818

1819

Run_Update(locEventLoop);

1820

1821

string locParticleNamesForHistX = "";

1822

for(size_t loc_i = 0; loc_i < dXPIDs.size(); ++loc_i)

1823

locParticleNamesForHistX += ParticleName_ROOT(dXPIDs[loc_i]);

1824

1825

string locParticleNamesForHistY = "";

1826

for(size_t loc_i = 0; loc_i < dYPIDs.size(); ++loc_i)

1827

locParticleNamesForHistY += ParticleName_ROOT(dYPIDs[loc_i]);

1828

1829

string locMassString = " Invariant Mass (GeV/c^{2});";

1830

1831

//CREATE THE HISTOGRAMS

1832

//Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock

1833

japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!

1834

{

1835

CreateAndChangeTo_ActionDirectory();

1836

1837

locHistName = "2DInvariantMass";

1838

locHistTitle = string(";") + locParticleNamesForHistX + locMassString + locParticleNamesForHistY + locMassString;

1839

dHist_2DInvaraintMass = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumXBins, dMinX, dMaxX, dNumYBins, dMinY, dMaxY);

1840

1841

//Return to the base directory

1842

ChangeTo_BaseDirectory();

1843

}

1844

japp->RootUnLock(); //RELEASE ROOT LOCK!!

1845

}

1846

1847

void DHistogramAction_2DInvariantMass::Run_Update(JEventLoop* locEventLoop)

1848

{

1849

DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());

1850

DGeometry *locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());

1851

locGeometry->GetTargetZ(dTargetZCenter);

1852

1853

locEventLoop->GetSingle(dAnalysisUtilities);

1854

1855

//BEAM BUNCH PERIOD

1856

vector<double> locBeamPeriodVector;

1857

locEventLoop->GetCalib("PHOTON_BEAM/RF/beam_period", locBeamPeriodVector);

1858

dBeamBunchPeriod = locBeamPeriodVector[0];

1859

}

1860

1861

bool DHistogramAction_2DInvariantMass::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)

1862

{

1863

vector<pair<double, double> > locMassesToFill;

1864

for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i)

1865

{

1866

const DReactionStep* locReactionStep = Get_Reaction()->Get_ReactionStep(loc_i);

1867

if((dStepIndex != -1) && (int(loc_i) != dStepIndex))

1868

continue;

1869

1870

//build all possible combinations of the included pids

1871

set<set<size_t> > locXIndexCombos = dAnalysisUtilities->Build_IndexCombos(locReactionStep, dXPIDs);

1872

set<set<size_t> > locYIndexCombos = dAnalysisUtilities->Build_IndexCombos(locReactionStep, dYPIDs);

1873

1874

//(double) loop over them

1875

set<set<size_t> >::iterator locXComboIterator = locXIndexCombos.begin();

1876

for(; locXComboIterator != locXIndexCombos.end(); ++locXComboIterator)

1877

{

1878

set<pair<const JObject*, unsigned int> > locXSourceObjects;

1879

DLorentzVector locXP4 = dAnalysisUtilities->Calc_FinalStateP4(Get_Reaction(), locParticleCombo, loc_i, *locXComboIterator, locXSourceObjects, Get_UseKinFitResultsFlag());

1880

1881

set<set<size_t> >::iterator locYComboIterator = locYIndexCombos.begin();

1882

for(; locYComboIterator != locYIndexCombos.end(); ++locYComboIterator)

1883

{

1884

set<pair<const JObject*, unsigned int> > locYSourceObjects;

1885

DLorentzVector locYP4 = dAnalysisUtilities->Calc_FinalStateP4(Get_Reaction(), locParticleCombo, loc_i, *locYComboIterator, locYSourceObjects, Get_UseKinFitResultsFlag());

1886

1887

if(locXSourceObjects == locYSourceObjects)

1888

continue; //the same!

1889

1890

set<set<pair<const JObject*, unsigned int> > > locAllSourceObjects;

1891

locAllSourceObjects.insert(locXSourceObjects);

1892

locAllSourceObjects.insert(locYSourceObjects);

1893

if(dPreviousSourceObjects.find(locAllSourceObjects) != dPreviousSourceObjects.end())

1894

continue; //dupe: already histed! (also, could be that the X/Y swapped combo has already been histed: don't double-count!

1895

dPreviousSourceObjects.insert(locAllSourceObjects);

1896

1897

locMassesToFill.push_back(pair<double, double>(locXP4.M(), locYP4.M()));

1898

}

1899

}

1900

}

1901

1902

// calculate accidental scaling factor (if needed)

1903

double locWeight = 1.;

1904

if(dSubtractAccidentals) {

1905

auto locIsFirstStepBeam = DAnalysis::Get_IsFirstStepBeam(Get_Reaction());

1906

if(locIsFirstStepBeam) {

1907

const DEventRFBunch* locEventRFBunch = locParticleCombo->Get_EventRFBunch();

1908

const DKinematicData* locBeamParticle = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle();

1909

1910

double locDeltaTRF = locBeamParticle->time() - (locEventRFBunch->dTime + (locBeamParticle->z() - dTargetZCenter)/29.9792458);

1911

1912

if(fabs(locDeltaTRF) > dBeamBunchPeriod/2.)

1913

locWeight = -1./(2.*Get_Reaction()->Get_NumPlusMinusRFBunches());

1914

}

1915

}

1916

1917

//FILL HISTOGRAMS

1918

//Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock

1919

//Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex

1920

Lock_Action();

1921

{

1922

if(dSubtractAccidentals) {

1923

for(size_t loc_i = 0; loc_i < locMassesToFill.size(); ++loc_i)

1924

dHist_2DInvaraintMass->Fill(locMassesToFill[loc_i].first, locMassesToFill[loc_i].second, locWeight);

1925

} else {

1926

for(size_t loc_i = 0; loc_i < locMassesToFill.size(); ++loc_i)

1927

dHist_2DInvaraintMass->Fill(locMassesToFill[loc_i].first, locMassesToFill[loc_i].second);

1928

}

1929

}

1930

Unlock_Action();

1931

1932

return true;

1933

}

1934

1935

void DHistogramAction_Dalitz::Initialize(JEventLoop* locEventLoop)

1936

{

1937

string locHistName, locHistTitle;

1938

1939

Run_Update(locEventLoop);

1940

1941

string locParticleNamesForHistX = "";

1942

for(size_t loc_i = 0; loc_i < dXPIDs.size(); ++loc_i)

1943

locParticleNamesForHistX += ParticleName_ROOT(dXPIDs[loc_i]);

1944

1945

string locParticleNamesForHistY = "";

1946

for(size_t loc_i = 0; loc_i < dYPIDs.size(); ++loc_i)

1947

locParticleNamesForHistY += ParticleName_ROOT(dYPIDs[loc_i]);

1948

1949

string locMassString = " Invariant Mass Squared (GeV/c^{2})^{2};";

1950

1951

//CREATE THE HISTOGRAMS

1952

//Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock

1953

japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!

1954

{

1955

CreateAndChangeTo_ActionDirectory();

1956

1957

locHistName = "DalitzPlot";

1958

locHistTitle = string(";") + locParticleNamesForHistX + locMassString + locParticleNamesForHistY + locMassString;

1959

dHist_DalitzPlot = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNumXBins, dMinX, dMaxX, dNumYBins, dMinY, dMaxY);

1960

1961

//Return to the base directory

1962

ChangeTo_BaseDirectory();

1963

}

1964

japp->RootUnLock(); //RELEASE ROOT LOCK!!

1965

}

1966

1967

void DHistogramAction_Dalitz::Run_Update(JEventLoop* locEventLoop)

1968

{

1969

locEventLoop->GetSingle(dAnalysisUtilities);

1970

}

1971

1972

bool DHistogramAction_Dalitz::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)

1973

{

1974

vector<pair<double, double> > locMassesToFill;

1975

for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i)

1976

{

1977

const DReactionStep* locReactionStep = Get_Reaction()->Get_ReactionStep(loc_i);

1978

if((dStepIndex != -1) && (int(loc_i) != dStepIndex))

1979

continue;

1980

1981

//build all possible combinations of the included pids

1982

set<set<size_t> > locXIndexCombos = dAnalysisUtilities->Build_IndexCombos(locReactionStep, dXPIDs);

1983

set<set<size_t> > locYIndexCombos = dAnalysisUtilities->Build_IndexCombos(locReactionStep, dYPIDs);

1984

1985

//(double) loop over them

1986

set<set<size_t> >::iterator locXComboIterator = locXIndexCombos.begin();

1987

for(; locXComboIterator != locXIndexCombos.end(); ++locXComboIterator)

1988

{

1989

set<pair<const JObject*, unsigned int> > locXSourceObjects;

1990

DLorentzVector locXP4 = dAnalysisUtilities->Calc_FinalStateP4(Get_Reaction(), locParticleCombo, loc_i, *locXComboIterator, locXSourceObjects, Get_UseKinFitResultsFlag());

1991

1992

set<set<size_t> >::iterator locYComboIterator = locYIndexCombos.begin();

1993

for(; locYComboIterator != locYIndexCombos.end(); ++locYComboIterator)

1994

{

1995

set<pair<const JObject*, unsigned int> > locYSourceObjects;

1996

DLorentzVector locYP4 = dAnalysisUtilities->Calc_FinalStateP4(Get_Reaction(), locParticleCombo, loc_i, *locYComboIterator, locYSourceObjects, Get_UseKinFitResultsFlag());

1997

1998

if(locXSourceObjects == locYSourceObjects)

1999

continue; //the same!

2000

2001

set<set<pair<const JObject*, unsigned int> > > locAllSourceObjects;

2002

locAllSourceObjects.insert(locXSourceObjects);

2003

locAllSourceObjects.insert(locYSourceObjects);

2004

if(dPreviousSourceObjects.find(locAllSourceObjects) != dPreviousSourceObjects.end())

2005

continue; //dupe: already histed! (also, could be that the X/Y swapped combo has already been histed: don't double-count!

2006

dPreviousSourceObjects.insert(locAllSourceObjects);

2007

2008

locMassesToFill.push_back(pair<double, double>(locXP4.M2(), locYP4.M2()));

2009

}

2010

}

2011

}

2012

2013

//FILL HISTOGRAMS

2014

//Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock

2015

//Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex

2016

Lock_Action();

2017

{

2018

for(size_t loc_i = 0; loc_i < locMassesToFill.size(); ++loc_i)

2019

dHist_DalitzPlot->Fill(locMassesToFill[loc_i].first, locMassesToFill[loc_i].second);

2020

}

2021

Unlock_Action();

2022

2023

return true;

2024

}

2025

2026

void DHistogramAction_KinFitResults::Initialize(JEventLoop* locEventLoop)

2027

{

2028

if (gPARMS->Exists("KINFIT:DEPENDENCE_HISTS")){

2029

bool locHistDependenceFlag = false;

2030

gPARMS->SetDefaultParameter("KINFIT:DEPENDENCE_HISTS", locHistDependenceFlag);

2031

gPARMS->GetParameter("KINFIT:DEPENDENCE_HISTS", locHistDependenceFlag);

2032

dHistDependenceFlag = locHistDependenceFlag;

2033

}

2034

2035

auto locReaction = Get_Reaction();

2036

DKinFitType locKinFitType = locReaction->Get_KinFitType();

2037

if(locKinFitType == d_NoFit)

2038

return;

2039

2040

//Get the DReactionVertexInfo for this reaction

2041

vector<const DReactionVertexInfo*> locReactionVertexInfos;

2042

locEventLoop->Get(locReactionVertexInfos);

2043

const DReactionVertexInfo* locReactionVertexInfo = nullptr;

2044

for(auto& locVertexInfo : locReactionVertexInfos)

2045

{

2046

auto locReactions = locVertexInfo->Get_Reactions();

2047

std::sort(locReactions.begin(), locReactions.end());

2048

if(!std::binary_search(locReactions.begin(), locReactions.end(), Get_Reaction()))

2049

continue;

2050

locReactionVertexInfo = locVertexInfo;

2051

break;

2052

}

2053

2054

dKinFitUtils = new DKinFitUtils_GlueX(locEventLoop);

2055

Run_Update(locEventLoop);

2056

2057

size_t locNumConstraints = 0, locNumUnknowns = 0;

2058

string locConstraintString = dKinFitUtils->Get_ConstraintInfo(locReactionVertexInfo, Get_Reaction(), locNumConstraints, locNumUnknowns);

2059

2060

size_t locNDF = locNumConstraints - locNumUnknowns;

2061

bool locIncludeBeamlineInVertexFitFlag = dKinFitUtils->Get_IncludeBeamlineInVertexFitFlag();

2062

2063

bool locP4IsFit = ((locKinFitType != d_VertexFit) && (locKinFitType != d_SpacetimeFit));

2064

bool locSpactimeIsFitFlag = (locKinFitType == d_SpacetimeFit) || (locKinFitType == d_P4AndSpacetimeFit);

2065

bool locVertexIsFitFlag = (locSpactimeIsFitFlag || (locKinFitType == d_VertexFit) || (locKinFitType == d_P4AndVertexFit));

2066

2067

//bool locIsInclusiveChannelFlag = Get_Reaction()->Get_IsInclusiveChannelFlag();

2068

//Below, should in theory check on whether to create pxyz pull histograms in the inclusive channel case

2069

//But, this is tricky: can have inclusive (no p4) but still have mass constraints (p4)

2070

2071

//CREATE THE HISTOGRAMS

2072

//Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock

2073

japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!

2074

{

2075

CreateAndChangeTo_ActionDirectory();

2076

2077

// Confidence Level

2078

string locHistName = "ConfidenceLevel";

2079

ostringstream locHistTitle;

2080

locHistTitle << "Kinematic Fit Constraints: " << locConstraintString << ";Confidence Level (" << locNumConstraints;

2081

locHistTitle << " Constraints, " << locNumUnknowns << " Unknowns: " << locNDF << "-C Fit);# Combos";

2082

dHist_ConfidenceLevel = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle.str(), dNumConfidenceLevelBins, 0.0, 1.0);

2083

2084

// Reduced chi^2 (chi^2/# d.o.f.)

2085

locHistName = "ChiSq";

2086

locHistTitle.str(""); // clear the ostringstream

2087

locHistTitle << "Kinematic Fit Constraints: " << locConstraintString << "; #chi^{2}/# d.o.f. (" << locNumConstraints;

2088

locHistTitle << " Constraints, " << locNumUnknowns << " Unknowns: " << locNDF << "-C Fit);# Combos";

2089

dHist_ChiSq = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle.str(), dNumChiSqBins, 0.0, dMaxChiSq);

2090

2091

//beam (pulls & conlev)

2092

Particle_t locInitialPID = Get_Reaction()->Get_ReactionStep(0)->Get_InitialPID();

2093

bool locBeamFlag = Get_IsFirstStepBeam(Get_Reaction());

2094

if(locBeamFlag)

2095

{

2096

auto locStepVertexInfo = locReactionVertexInfo->Get_StepVertexInfo(0);

2097

auto locFullConstrainParticles = locStepVertexInfo->Get_FullConstrainParticles(true, d_InitialState, d_AllCharges, false);

2098

bool locIsInVertexFitFlag = locVertexIsFitFlag && locIncludeBeamlineInVertexFitFlag && !locFullConstrainParticles.empty();

2099

bool locIsChargedFlag = (ParticleCharge(locInitialPID) != 0);

2100

2101

if(locP4IsFit || locIsInVertexFitFlag)

2102

{

2103

string locFullROOTName = string("Beam ") + ParticleName_ROOT(locInitialPID);

2104

2105

if(dHistDependenceFlag)

2106

{

2107

//Conlev

2108

locHistName = "ConfidenceLevel_VsBeamEnergy";

2109

locHistTitle.str("");

2110

locHistTitle << "Kinematic Fit Constraints: " << locConstraintString << ";Beam Energy (GeV);Confidence Level (" << locNumConstraints;

2111

locHistTitle << " Constraints, " << locNumUnknowns << " Unknowns: " << locNDF << "-C Fit)";

2112

dHistMap_ConfidenceLevel_VsP[pair<int, Particle_t>(-1, locInitialPID)] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle.str(), dNum2DBeamEBins, dMinBeamE, dMaxBeamE, dNum2DConfidenceLevelBins, 0.0, 1.0);

2113

}

2114

2115

//Pulls

2116

CreateAndChangeTo_Directory("Beam", "Beam");

2117

map<DKinFitPullType, TH1I*> locParticlePulls;

2118

map<DKinFitPullType, TH2I*> locParticlePullsVsP, locParticlePullsVsTheta, locParticlePullsVsPhi;

2119

Create_ParticlePulls(locFullROOTName, locIsChargedFlag, locIsInVertexFitFlag, false, -1, locInitialPID);

2120

2121

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

2122

}

2123

}

2124

2125

//final particles

2126

for(size_t loc_i = 0; loc_i < Get_Reaction()->Get_NumReactionSteps(); ++loc_i)

2127

{

2128

const DReactionStep* locReactionStep = Get_Reaction()->Get_ReactionStep(loc_i);

2129

auto locStepVertexInfo = locReactionVertexInfo->Get_StepVertexInfo(loc_i);

2130

auto locFullConstrainParticles = locStepVertexInfo->Get_FullConstrainParticles(true, d_FinalState, d_AllCharges, false);

2131

ostringstream locStepName;

2132

locStepName << "Step" << loc_i << "__" << DAnalysis::Get_StepName(locReactionStep, true, false);

2133

string locStepROOTName = DAnalysis::Get_StepName(locReactionStep, true, true);

2134

2135

//Conlev

2136

auto locPIDs = Get_Reaction()->Get_FinalPIDs(loc_i, false, false, d_AllCharges, false);

2137

for(auto locPID : locPIDs)

2138

{

2139

bool locIsInVertexFitFlag = locVertexIsFitFlag && locStepVertexInfo->Get_FittableVertexFlag();

2140

2141

bool locIsNeutralShowerFlag = (locIsInVertexFitFlag && (ParticleCharge(locPID) == 0));

2142

if((ParticleMass(locPID) > 0.0) && !locSpactimeIsFitFlag)

2143

locIsNeutralShowerFlag = false; //massive shower momentum is defined by t, which isn't fit: use particle

2144

if(!locP4IsFit && !locIsInVertexFitFlag)

2145

continue; //p4 is not fit, and this is not in a vertex fit: no pulls

2146

if(locIsNeutralShowerFlag && !locP4IsFit)

2147

continue; //vertex-only fit: neutral shower does not constrain: no pulls

2148

2149

string locParticleName = ParticleType(locPID);

2150

string locParticleROOTName = ParticleName_ROOT(locPID);

2151

2152

if(dHistDependenceFlag)

2153

{

2154

//vs p

2155

locHistName = string("ConfidenceLevel_Vs") + locParticleName + string("P");

2156

locHistTitle.str("");

2157

locHistTitle << "Kinematic Fit Constraints: " << locConstraintString << ";" << locParticleROOTName << " p (GeV/c);Confidence Level (" << locNumConstraints;

2158

locHistTitle << " Constraints, " << locNumUnknowns << " Unknowns: " << locNDF << "-C Fit)";

2159

dHistMap_ConfidenceLevel_VsP[pair<int, Particle_t>(loc_i, locPID)] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle.str(), dNum2DPBins, dMinP, dMaxP, dNum2DConfidenceLevelBins, 0.0, 1.0);

2160

2161

//vs theta

2162

locHistName = string("ConfidenceLevel_Vs") + locParticleName + string("Theta");

2163

locHistTitle.str("");

2164

locHistTitle << "Kinematic Fit Constraints: " << locConstraintString << ";" << locParticleROOTName << " #theta#circ;Confidence Level (" << locNumConstraints;

2165

locHistTitle << " Constraints, " << locNumUnknowns << " Unknowns: " << locNDF << "-C Fit)";

2166

dHistMap_ConfidenceLevel_VsTheta[pair<int, Particle_t>(loc_i, locPID)] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle.str(), dNum2DThetaBins, dMinTheta, dMaxTheta, dNum2DConfidenceLevelBins, 0.0, 1.0);

2167

2168

//vs phi

2169

locHistName = string("ConfidenceLevel_Vs") + locParticleName + string("Phi");

2170

locHistTitle.str("");

2171

locHistTitle << "Kinematic Fit Constraints: " << locConstraintString << ";" << locParticleROOTName << " #phi#circ;Confidence Level (" << locNumConstraints;

2172

locHistTitle << " Constraints, " << locNumUnknowns << " Unknowns: " << locNDF << "-C Fit)";

2173

dHistMap_ConfidenceLevel_VsPhi[pair<int, Particle_t>(loc_i, locPID)] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle.str(), dNum2DPhiBins, dMinPhi, dMaxPhi, dNum2DConfidenceLevelBins, 0.0, 1.0);

2174

}

2175

}

2176

2177

//pulls

2178

bool locFilledFlag = false;

2179

for(auto locPID : locPIDs)

2180

{

2181

bool locIsInVertexFitFlag = locVertexIsFitFlag && locStepVertexInfo->Get_FittableVertexFlag();

2182

bool locIsChargedFlag = (ParticleCharge(locPID) != 0);

2183

2184

bool locIsNeutralShowerFlag = (locIsInVertexFitFlag && (ParticleCharge(locPID) == 0));

2185

if((ParticleMass(locPID) > 0.0) && !locSpactimeIsFitFlag)

2186

locIsNeutralShowerFlag = false; //massive shower momentum is defined by t, which isn't fit: use particle

2187

if(!locP4IsFit && !locIsInVertexFitFlag)

2188

continue; //p4 is not fit, and this is not in a vertex fit: no pulls

2189

if(locIsNeutralShowerFlag && !locP4IsFit)

2190

continue; //vertex-only fit: neutral shower does not constrain: no pulls

2191

2192

if(!locFilledFlag) //first call

2193

{

2194

locFilledFlag = true;

2195

CreateAndChangeTo_Directory(locStepName.str(), locStepName.str());

2196

}

2197

2198

string locParticleName = ParticleType(locPID);

2199

CreateAndChangeTo_Directory(locParticleName, locParticleName);

2200

2201

string locParticleROOTName = ParticleName_ROOT(locPID);

2202

string locFullROOTName = locParticleROOTName + string(", ") + locStepROOTName;

2203

Create_ParticlePulls(locFullROOTName, locIsChargedFlag, locIsInVertexFitFlag, locIsNeutralShowerFlag, loc_i, locPID);

2204

2205

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

2206

} //end of particle loop

2207

2208

if(locFilledFlag)

2209

gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->cd("..");

2210

} //end of step loop

2211

2212

//Return to the base directory

2213

ChangeTo_BaseDirectory();

2214

}

2215

japp->RootUnLock(); //RELEASE ROOT LOCK!!

2216

}

2217

2218

void DHistogramAction_KinFitResults::Run_Update(JEventLoop* locEventLoop)

2219

{

2220

locEventLoop->GetSingle(dAnalysisUtilities);

2221

dKinFitUtils->Set_RunDependent_Data(locEventLoop);

2222

}

2223

2224

void DHistogramAction_KinFitResults::Create_ParticlePulls(string locFullROOTName, bool locIsChargedFlag, bool locIsInVertexFitFlag, bool locIsNeutralShowerFlag, int locStepIndex, Particle_t locPID)

2225

{

2226

auto locParticlePair = std::make_pair(locStepIndex, locPID);

2227

2228

DKinFitType locKinFitType = Get_Reaction()->Get_KinFitType();

2229

bool locP4IsFit = ((locKinFitType != d_VertexFit) && (locKinFitType != d_SpacetimeFit));

Assuming 'locKinFitType' is equal to d_VertexFit

→

2230

2231

//Determine pull types

2232

set<pair<DKinFitPullType, pair<string, string> > > locPullTypes;

2233

2234

//p4 pulls:

2235

string locHistName, locHistTitle;

2236

if(locIsNeutralShowerFlag && locP4IsFit) //neutral shower not in a p4-only fit

2237

locPullTypes.insert(pair<DKinFitPullType, pair<string, string> >(d_EPull, pair<string, string>("Pull_E", "E Pull")));

2238

else if(!locIsNeutralShowerFlag && (locP4IsFit || locIsInVertexFitFlag))

←

Taking false branch

→

2239

{

2240

//all detected particles (except neutral showers) have p3 used in p4 fits and vertex fits

2241

//however, don't include if the particles aren't actually in the fit (vertex-only, and too few particles at that vertex to constrain)

2242

locPullTypes.insert(pair<DKinFitPullType, pair<string, string> >(d_PxPull, pair<string, string>("Pull_Px", "p_{x} Pull")));

2243

locPullTypes.insert(pair<DKinFitPullType, pair<string, string> >(d_PyPull, pair<string, string>("Pull_Py", "p_{y} Pull")));

2244

locPullTypes.insert(pair<DKinFitPullType, pair<string, string> >(d_PzPull, pair<string, string>("Pull_Pz", "p_{z} Pull")));

2245

}

2246

2247

//vertex pulls:

2248

if((locIsNeutralShowerFlag && locP4IsFit) || (locIsChargedFlag && locIsInVertexFitFlag))

2249

{

2250

locPullTypes.insert(pair<DKinFitPullType, pair<string, string> >(d_XxPull, pair<string, string>("Pull_Xx", "x_{x} Pull")));

2251

locPullTypes.insert(pair<DKinFitPullType, pair<string, string> >(d_XyPull, pair<string, string>("Pull_Xy", "x_{y} Pull")));

2252

locPullTypes.insert(pair<DKinFitPullType, pair<string, string> >(d_XzPull, pair<string, string>("Pull_Xz", "x_{z} Pull")));

2253

}

2254

2255

//time pulls:

2256

if(locIsInVertexFitFlag && ((locKinFitType == d_SpacetimeFit) || (locKinFitType == d_P4AndSpacetimeFit)))

2257

locPullTypes.insert(pair<DKinFitPullType, pair<string, string> >(d_TPull, pair<string, string>("Pull_T", "t Pull")));

2258

2259

bool locIsBeamFlag = (locFullROOTName.substr(0, 4) == "Beam");

2260

int locNum2DPBins = locIsBeamFlag ? dNum2DBeamEBins : dNum2DPBins;

←

'?' condition is false

→

2261

double locMinP = locIsBeamFlag ? dMinBeamE : dMinP;

←

'?' condition is false

→

2262

double locMaxP = locIsBeamFlag ? dMaxBeamE : dMaxP;

←

'?' condition is false

→

2263

2264

//Create histograms

2265

set<pair<DKinFitPullType, pair<string, string> > >::iterator locIterator = locPullTypes.begin();

2266

for(; locIterator != locPullTypes.end(); ++locIterator)

←

Loop condition is true. Entering loop body

→

2267

{

2268

pair<string, string> locStrings = (*locIterator).second;

2269

auto locPullType = (*locIterator).first;

2270

2271

//1D

2272

locHistName = locStrings.first;

2273

locHistTitle = locFullROOTName + string(";") + locStrings.second + string(";# Combos");

2274

dHistMap_Pulls[locParticlePair][locPullType] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPullBins, dMinPull, dMaxPull);

2275

2276

int locNumDeltaBins;

2277

double locMaxDeltaRange;

←

'locMaxDeltaRange' declared without an initial value

→

2278

Get_DeltaBinningParams(locPullType, false, locNumDeltaBins, locMaxDeltaRange);

←

Calling 'DHistogramAction_KinFitResults::Get_DeltaBinningParams'

→

←

Returning from 'DHistogramAction_KinFitResults::Get_DeltaBinningParams'

→

2279

Get_HistNameTitle(locPullType, locFullROOTName, 0, locHistName, locHistTitle);

2280

dHistMap_Deltas[locParticlePair][locPullType] = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, locNumDeltaBins, -1.0*locMaxDeltaRange, locMaxDeltaRange);

←

The right operand of '*' is a garbage value

2281

2282

if(!dHistDependenceFlag)

2283

continue;

2284

Get_DeltaBinningParams(locPullType, true, locNumDeltaBins, locMaxDeltaRange);

2285

2286

//vs p

2287

locHistName = locStrings.first + string("_VsP");

2288

locHistTitle = locFullROOTName + string(";p (GeV/c);") + locStrings.second;

2289

dHistMap_PullsVsP[locParticlePair][locPullType] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, locNum2DPBins, locMinP, locMaxP, dNumPullBins, dMinPull, dMaxPull);

2290

Get_HistNameTitle(locPullType, locFullROOTName, 1, locHistName, locHistTitle);

2291

dHistMap_DeltasVsP[locParticlePair][locPullType] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, locNum2DPBins, locMinP, locMaxP, locNumDeltaBins, -1.0*locMaxDeltaRange, locMaxDeltaRange);

2292

if(locIsBeamFlag)

2293

continue;

2294

2295

//vs theta

2296

locHistName = locStrings.first + string("_VsTheta");

2297

locHistTitle = locFullROOTName + string(";#theta#circ;") + locStrings.second;

2298

dHistMap_PullsVsTheta[locParticlePair][locPullType] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, dNumPullBins, dMinPull, dMaxPull);

2299

Get_HistNameTitle(locPullType, locFullROOTName, 2, locHistName, locHistTitle);

2300

dHistMap_DeltasVsTheta[locParticlePair][locPullType] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DThetaBins, dMinTheta, dMaxTheta, locNumDeltaBins, -1.0*locMaxDeltaRange, locMaxDeltaRange);

2301

2302

//vs phi

2303

locHistName = locStrings.first + string("_VsPhi");

2304

locHistTitle = locFullROOTName + string(";#phi#circ;") + locStrings.second;

2305

dHistMap_PullsVsPhi[locParticlePair][locPullType] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, dNumPullBins, dMinPull, dMaxPull);

2306

Get_HistNameTitle(locPullType, locFullROOTName, 3, locHistName, locHistTitle);

2307

dHistMap_DeltasVsPhi[locParticlePair][locPullType] = GetOrCreate_Histogram<TH2I>(locHistName, locHistTitle, dNum2DPhiBins, dMinPhi, dMaxPhi, locNumDeltaBins, -1.0*locMaxDeltaRange, locMaxDeltaRange);

2308

}

2309

}

2310

2311

void DHistogramAction_KinFitResults::Get_HistNameTitle(DKinFitPullType locPullType, string locFullROOTName, int locVsKey, string& locHistName, string& locHistTitle)

2312

{

2313

string locNameParam, locTitleParam;

2314

if(locPullType == d_EPull)

2315

{

2316

locNameParam = "ShowerE";

2317

locTitleParam = "E";

2318

}

2319

else if(locPullType == d_PxPull) //delta-theta!!

2320

{

2321

locNameParam = "Theta";

2322

locTitleParam = "#theta#circ";

2323

}

2324

else if(locPullType == d_PyPull) //delta-phi!!

2325

{

2326

locNameParam = "Phi";

2327

locTitleParam = "#phi#circ";

2328

}

2329

else if(locPullType == d_PzPull) //delta-p/p!!

2330

{

2331

locNameParam = "POverP";

2332

locTitleParam = "p/p";

2333

}

2334

else if(locPullType == d_XxPull)

2335

{

2336

locNameParam = "X";

2337

locTitleParam = "Vertex-X (cm)";

2338

}

2339

else if(locPullType == d_XyPull)

2340

{

2341

locNameParam = "Y";

2342

locTitleParam = "Vertex-Y (cm)";

2343

}

2344

else if(locPullType == d_XzPull)

2345

{

2346

locNameParam = "Z";

2347

locTitleParam = "Vertex-Z (cm)";

2348

}

2349

else if(locPullType == d_TPull)

2350

{

2351

locNameParam = "T";

2352

locTitleParam = "T (ns)";

2353

}

2354

2355

if(locVsKey == 0) //1D

2356

{

2357

locHistName = string("Delta") + locNameParam;

2358

locHistTitle = locFullROOTName + string(";#Delta") + locTitleParam;

2359

}

2360

else if(locVsKey == 1) //p

2361

{

2362

locHistName = string("Delta") + locNameParam + string("_VsP");

2363

locHistTitle = locFullROOTName + string(";p (GeV/c);#Delta") + locTitleParam;

2364

}

2365

else if(locVsKey == 2) //theta

2366

{

2367

locHistName = string("Delta") + locNameParam + string("_VsTheta");

2368

locHistTitle = locFullROOTName + string(";#theta#circ;#Delta") + locTitleParam;

2369

}

2370

else if(locVsKey == 3) //phi

2371

{

2372

locHistName = string("Delta") + locNameParam + string("_VsPhi");

2373

locHistTitle = locFullROOTName + string(";#phi#circ;#Delta") + locTitleParam;

2374

}

2375

}

2376

2377

void DHistogramAction_KinFitResults::Get_DeltaBinningParams(DKinFitPullType locPullType, bool loc2DFlag, int& locNumBins, double& locMax)

2378

{

2379

if(locPullType == d_EPull)

←

Assuming 'locPullType' is not equal to d_EPull

→

←

Taking false branch

→

2380

{

2381

locNumBins = loc2DFlag ? dNum2DDeltaShowerEBins : dNumDeltaShowerEBins;

2382

locMax = dMaxDeltaShowerE;

2383

}

2384

else if(locPullType == d_PxPull) //delta-theta!!

←

Assuming 'locPullType' is not equal to d_PxPull

→

←

Taking false branch

→

2385

{

2386

locNumBins = loc2DFlag ? dNum2DDeltaThetaBins : dNumDeltaThetaBins;

2387

locMax = dMaxDeltaTheta;

2388

}

2389

else if(locPullType == d_PyPull) //delta-phi!!

←

Assuming 'locPullType' is not equal to d_PyPull

→

←

Taking false branch

→

2390

{

2391

locNumBins = loc2DFlag ? dNum2DDeltaPhiBins : dNumDeltaPhiBins;

2392

locMax = dMaxDeltaPhi;

2393

}

2394

else if(locPullType == d_PzPull) //delta-p/p!!

←

Assuming 'locPullType' is not equal to d_PzPull

→

←

Taking false branch

→

2395

{

2396

locNumBins = loc2DFlag ? dNum2DDeltaPOverPBins : dNumDeltaPOverPBins;

2397

locMax = dMaxDeltaPOverP;

2398

}

2399

else if((locPullType == d_XxPull) || (locPullType == d_XyPull))

←

Assuming 'locPullType' is not equal to d_XxPull

→

←

Assuming 'locPullType' is not equal to d_XyPull

→

←

Taking false branch

→

2400

{

2401

locNumBins = loc2DFlag ? dNum2DDeltaVertexXYBins : dNumDeltaVertexXYBins;

2402

locMax = dMaxDeltaVertexXY;

2403

}

2404

else if(locPullType == d_XzPull)

←

Assuming 'locPullType' is not equal to d_XzPull

→

←

Taking false branch

→

2405

{

2406

locNumBins = loc2DFlag ? dNum2DDeltaVertexZBins : dNumDeltaVertexZBins;

2407

locMax = dMaxDeltaVertexZ;

2408

}

2409

else if(locPullType == d_TPull)

←

Assuming 'locPullType' is not equal to d_TPull

→

←

Taking false branch

→

2410

{

2411

locNumBins = loc2DFlag ? dNum2DDeltaTBins : dNumDeltaTBins;

2412

locMax = dMaxDeltaT;

2413

}

2414

}

2415

2416

bool DHistogramAction_KinFitResults::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)

2417

{

2418

//kinfit results are unique for each DParticleCombo: no need to check for duplicates

2419

const DKinFitResults* locKinFitResults = locParticleCombo->Get_KinFitResults();

2420

if(locKinFitResults == NULL__null)

2421

return true;

2422

2423

// Get Confidence Level

2424

double locConfidenceLevel = locKinFitResults->Get_ConfidenceLevel();

2425

2426

// Get reduced chi^2

2427

double locChiSq = locKinFitResults->Get_ChiSq()/locKinFitResults->Get_NDF() ;

2428

2429

// Get Pulls

2430

map<const JObject*, map<DKinFitPullType, double> > locPulls; //DKinematicData is the MEASURED particle

2431

locKinFitResults->Get_Pulls(locPulls);

2432

2433

Particle_t locInitPID = Get_Reaction()->Get_ReactionStep(0)->Get_InitialPID();

2434

bool locBeamFlag = DAnalysis::Get_IsFirstStepBeam(Get_Reaction());

2435

2436

//FILL HISTOGRAMS

2437

//Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock

2438

//Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex

2439

Lock_Action();

2440

{

2441

dHist_ConfidenceLevel->Fill(locConfidenceLevel);

2442

dHist_ChiSq->Fill(locChiSq);

2443

2444

// beam

2445

if(locBeamFlag)

2446

{

2447

const DKinematicData* locKinematicData = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle_Measured();

2448

map<DKinFitPullType, double> locParticlePulls = locPulls[locKinematicData];

2449

if(!locParticlePulls.empty())

2450

{

2451

pair<int, Particle_t> locParticlePair(-1, locInitPID);

2452

DVector3 locMomentum = locKinematicData->momentum();

2453

double locP = locMomentum.Mag();

2454

2455

//con lev

2456

if(dHistDependenceFlag)

2457

dHistMap_ConfidenceLevel_VsP[locParticlePair]->Fill(locP, locConfidenceLevel);

2458

2459

//pulls

2460

if(locConfidenceLevel >= dPullHistConfidenceLevelCut)

2461

{

2462

map<DKinFitPullType, double>::iterator locIterator = locParticlePulls.begin();

2463

for(; locIterator != locParticlePulls.end(); ++locIterator)

2464

{

2465

dHistMap_Pulls[locParticlePair][locIterator->first]->Fill(locIterator->second);

2466

if(dHistDependenceFlag)

2467

dHistMap_PullsVsP[locParticlePair][locIterator->first]->Fill(locP, locIterator->second);

2468

}

2469

}

2470

2471

//deltas

2472

double locDeltaPOverP = (locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle()->pmag() - locP)/locP;

2473

dHistMap_Deltas[locParticlePair][d_PzPull]->Fill(locDeltaPOverP);

2474

if(dHistDependenceFlag)

2475

dHistMap_DeltasVsP[locParticlePair][d_PzPull]->Fill(locP, locDeltaPOverP);

2476

}

2477

}

2478

2479

// final particle pulls

2480

for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i)

2481

{

2482

const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(loc_i);

2483

auto locParticles = locParticleComboStep->Get_FinalParticles_Measured(Get_Reaction()->Get_ReactionStep(loc_i));

2484

auto locKinFitParticles = locParticleComboStep->Get_FinalParticles(Get_Reaction()->Get_ReactionStep(loc_i), false, false);

2485

for(size_t loc_j = 0; loc_j < locParticles.size(); ++loc_j)

2486

{

2487

auto& locParticle = locParticles[loc_j];

2488

auto& locKinFitParticle = locKinFitParticles[loc_j];

2489

2490

//get pulls for this particle

2491

map<DKinFitPullType, double> locParticlePulls;

2492

map<const JObject*, map<DKinFitPullType, double> >::iterator locParticleIterator = locPulls.find(locParticle);

2493

if(locParticleIterator != locPulls.end())

2494

locParticlePulls = locParticleIterator->second;

2495

else //is neutral shower

2496

{

2497

auto locNeutralHypo = dynamic_cast<const DNeutralParticleHypothesis*>(locParticle);

2498

locParticlePulls = locPulls[locNeutralHypo->Get_NeutralShower()];

2499

}

2500

if(locParticlePulls.empty())

2501

continue;

2502

2503

pair<int, Particle_t> locParticlePair(loc_i, locParticle->PID());

2504

DVector3 locMomentum = locParticle->momentum();

2505

double locP = locMomentum.Mag();

2506

double locTheta = locMomentum.Theta()*180.0/TMath::Pi();

2507

double locPhi = locMomentum.Phi()*180.0/TMath::Pi();

2508

2509

//deltas //loop used to easily figure out what params changed

2510

for(auto locIterator = locParticlePulls.begin(); locIterator != locParticlePulls.end(); ++locIterator)

2511

{

2512

auto locPullType = locIterator->first;

2513

auto locDelta = Get_Delta(locPullType, locParticle, locKinFitParticle);

2514

dHistMap_Deltas[locParticlePair][locPullType]->Fill(locDelta);

2515

if(!dHistDependenceFlag)

2516

continue;

2517

2518

dHistMap_DeltasVsP[locParticlePair][locPullType]->Fill(locP, locDelta);

2519

dHistMap_DeltasVsTheta[locParticlePair][locPullType]->Fill(locTheta, locDelta);

2520

dHistMap_DeltasVsPhi[locParticlePair][locPullType]->Fill(locPhi, locDelta);

2521

}

2522

2523

//con lev

2524

if(dHistDependenceFlag)

2525

{

2526

dHistMap_ConfidenceLevel_VsP[locParticlePair]->Fill(locP, locConfidenceLevel);

2527

dHistMap_ConfidenceLevel_VsTheta[locParticlePair]->Fill(locTheta, locConfidenceLevel);

2528

dHistMap_ConfidenceLevel_VsPhi[locParticlePair]->Fill(locPhi, locConfidenceLevel);

2529

}

2530

if(locConfidenceLevel < dPullHistConfidenceLevelCut)

2531

continue;

2532

2533

//fill histograms

2534

map<DKinFitPullType, double>::iterator locIterator = locParticlePulls.begin();

2535

for(; locIterator != locParticlePulls.end(); ++locIterator)

2536

{

2537

dHistMap_Pulls[locParticlePair][locIterator->first]->Fill(locIterator->second);

2538

if(!dHistDependenceFlag)

2539

continue;

2540

dHistMap_PullsVsP[locParticlePair][locIterator->first]->Fill(locP, locIterator->second);

2541

dHistMap_PullsVsTheta[locParticlePair][locIterator->first]->Fill(locTheta, locIterator->second);

2542

dHistMap_PullsVsPhi[locParticlePair][locIterator->first]->Fill(locPhi, locIterator->second);

2543

}

2544

}

2545

}

2546

}

2547

Unlock_Action();

2548

2549

return true;

2550

}

2551

2552

double DHistogramAction_KinFitResults::Get_Delta(DKinFitPullType locPullType, const DKinematicData* locMeasured, const DKinematicData* locKinFit)

2553

{

2554

//kinfit - measured

2555

if(locPullType == d_EPull)

2556

return locKinFit->energy() - locMeasured->energy();

2557

else if(locPullType == d_PxPull) //delta-theta!!

2558

return (locKinFit->momentum().Theta() - locMeasured->momentum().Theta())*180.0/TMath::Pi();

2559

else if(locPullType == d_PyPull) //delta-phi!!

2560

{

2561

auto locDeltaPhi = (locKinFit->momentum().Phi() - locMeasured->momentum().Phi())*180.0/TMath::Pi();

2562

while(locDeltaPhi > 180.0)

2563

locDeltaPhi -= 360.0;

2564

while(locDeltaPhi < -180.0)

2565

locDeltaPhi += 360.0;

2566

return locDeltaPhi;

2567

}

2568

else if(locPullType == d_PzPull) //delta-p-over-p!!

2569

{

2570

auto locP = locMeasured->momentum().Mag();

2571

return (locKinFit->momentum().Mag() - locP)/locP;

2572

}

2573

else if(locPullType == d_XxPull)

2574

return locKinFit->position().X() - locMeasured->position().X();

2575

else if(locPullType == d_XyPull)

2576

return locKinFit->position().Y() - locMeasured->position().Y();

2577

else if(locPullType == d_XzPull)

2578

return locKinFit->position().Z() - locMeasured->position().Z();

2579

else if(locPullType == d_TPull)

2580

return locKinFit->time() - locMeasured->time();

2581

else

2582

return 0.0;

2583

}

2584

2585

void DHistogramAction_MissingTransverseMomentum::Initialize(JEventLoop* locEventLoop)

2586

{

2587

string locHistName, locHistTitle;

2588

double locPtPerBin = 1000.0*(dMaxPt - dMinPt)/((double)dNumPtBins);

2589

2590

Run_Update(locEventLoop);

2591

2592

//CREATE THE HISTOGRAMS

2593

//Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock

2594

japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!

2595

{

2596

CreateAndChangeTo_ActionDirectory();

2597

2598

locHistName = "MissingTransverseMomentum";

2599

ostringstream locStream;

2600

locStream << locPtPerBin;

2601

locHistTitle = string(";") + string(" Missing Transverse Momentum (GeV/c);# Combos / ") + locStream.str() + string(" MeV/c");

2602

dHist_MissingTransverseMomentum = GetOrCreate_Histogram<TH1I>(locHistName, locHistTitle, dNumPtBins, dMinPt, dMaxPt);

2603

2604

//Return to the base directory

2605

ChangeTo_BaseDirectory();

2606

}

2607

japp->RootUnLock(); //RELEASE ROOT LOCK!!

2608

}

2609

2610

void DHistogramAction_MissingTransverseMomentum::Run_Update(JEventLoop* locEventLoop)

2611

{

2612

vector<const DAnalysisUtilities*> locAnalysisUtilitiesVector;

2613

locEventLoop->Get(locAnalysisUtilitiesVector);

2614

dAnalysisUtilities = locAnalysisUtilitiesVector[0];

2615

}

2616

2617

bool DHistogramAction_MissingTransverseMomentum::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)

2618

{

2619

set<pair<const JObject*, unsigned int> > locSourceObjects;

2620

DLorentzVector locFinalStateP4 = dAnalysisUtilities->Calc_FinalStateP4(Get_Reaction(), locParticleCombo, 0, locSourceObjects, Get_UseKinFitResultsFlag()); // Use step '0'

2621

2622

if(dPreviousSourceObjects.find(locSourceObjects) != dPreviousSourceObjects.end())

2623

return true; //dupe: already histed!

2624

dPreviousSourceObjects.insert(locSourceObjects);

2625

2626

double locMissingTransverseMomentum = locFinalStateP4.Pt();

2627

2628

//FILL HISTOGRAMS

2629

//Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock

2630

//Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex

2631

Lock_Action();

2632

{

2633

dHist_MissingTransverseMomentum->Fill(locMissingTransverseMomentum);

2634

}

2635

Unlock_Action();

2636

2637

return true;

2638

}

2639