libraries/ANALYSIS/DKinFitResults

Bug Summary

File:	libraries/ANALYSIS/DKinFitResults_factory.cc
Location:	line 411, column 58
Description:	Called C++ object pointer is null

Annotated Source Code

#ifdef VTRACE

#include "vt_user.h"

#endif

#include "DKinFitResults_factory.h"

//------------------

// init

//------------------

jerror_t DKinFitResults_factory::init(void)

{

dDebugLevel = 0;

dKinFitDebugLevel = 0;

dLinkVerticesFlag = true;

return NOERROR;

}

//------------------

// brun

//------------------

jerror_t DKinFitResults_factory::brun(jana::JEventLoop* locEventLoop, int runnumber)

{

vector<const DAnalysisUtilities*> locAnalysisUtilitiesVector;

locEventLoop->Get(locAnalysisUtilitiesVector);

if(locAnalysisUtilitiesVector.empty())

{

_DBG_std::cerr<<"libraries/ANALYSIS/DKinFitResults_factory.cc"
<<":"<<27<<" "<<"Unable to get a DAnalysisUtilities object!"<<endl;

return RESOURCE_UNAVAILABLE;

}

dAnalysisUtilities = locAnalysisUtilitiesVector[0];

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

const DMagneticFieldMap* locMagneticFieldMap = locApplication->GetBfield(runnumber);

dTargetZCenter = 65.0;

// Get Target parameters from XML

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

if(locGeometry != NULL__null)

locGeometry->GetTargetZ(dTargetZCenter);

double locBx, locBy, locBz;

locMagneticFieldMap->GetField(0.0, 0.0, dTargetZCenter, locBx, locBy, locBz);

TVector3 locBField(locBx, locBy, locBz);

if(locBField.Mag() > 0.0)

dKinFitter.Set_BField(locMagneticFieldMap);

gPARMS->SetDefaultParameter("KINFIT:KINFITDEBUGLEVEL", dKinFitDebugLevel);

gPARMS->SetDefaultParameter("KINFIT:DEBUGLEVEL", dDebugLevel);

gPARMS->SetDefaultParameter("KINFIT:LINKVERTICES", dLinkVerticesFlag);

dKinFitter.Set_DebugLevel(dKinFitDebugLevel);

dKinFitter.Set_LinkVerticesFlag(dLinkVerticesFlag);

// Get # of DReactions:

// Get list of factories and find all the ones producing

// DReaction objects. (A simpler way to do this would be to

// just use locEventLoop->Get(...), but then only one plugin could

// be used at a time.)

size_t locNumReactions = 0;

vector<JFactory_base*> locFactories = locEventLoop->GetFactories();

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

{

JFactory<DReaction>* locFactory = dynamic_cast<JFactory<DReaction>* >(locFactories[loc_i]);

if(locFactory == NULL__null)

continue;

if(string(locFactory->Tag()) == "Thrown")

continue;

// Found a factory producing DReactions. The reaction objects are

// produced at the init stage and are persistent through all event

// processing so we can grab the list here and append it to our

// overall list.

vector<const DReaction*> locReactionsSubset;

locFactory->Get(locReactionsSubset);

locNumReactions += locReactionsSubset.size();

}

//set pool sizes

unsigned int locExpectedNumCombos = 50;

dKinFitter.Set_MaxKinFitParticlePoolSize(locNumReactions*locExpectedNumCombos*6);

dKinFitter.Set_MaxKinFitConstraintVertexPoolSize(locNumReactions*locExpectedNumCombos*3);

dKinFitter.Set_MaxKinFitConstraintSpacetimePoolSize(locNumReactions*locExpectedNumCombos*3);

dKinFitter.Set_MaxKinFitConstraintP4PoolSize(locNumReactions*locExpectedNumCombos*3);

dKinFitter.Set_MaxMatrixDSymPoolSize(locNumReactions*locExpectedNumCombos*6);

dKinFitter.Set_MaxLargeMatrixDSymPoolSize(locNumReactions*locExpectedNumCombos);

dKinFitter.Preallocate_MatrixMemory();

return NOERROR;

}

void DKinFitResults_factory::Reset_NewEvent(void)

{

dKinFitter.Reset_NewEvent();

}

//------------------

// evnt

//------------------

jerror_t DKinFitResults_factory::evnt(JEventLoop* locEventLoop, int eventnumber)

100

{

101

#ifdef VTRACE

102

VT_TRACER("DKinFitResults_factory::evnt()");

103

#endif

104

105

dPreviouslyFailedFits.clear();

106

107

//perform all of the analysis steps that don't need the kinematic fit results (saves time by reducing #kinfits)

108

vector<const DAnalysisResults*> locAnalysisResultsVector;

109

locEventLoop->Get(locAnalysisResultsVector, "PreKinFit");

110

111

//get all of the ParticleCombos that survive the cuts

112

deque<const DParticleCombo*> locParticleCombos, locSurvivedParticleCombos;

113

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

114

{

115

locAnalysisResultsVector[loc_i]->Get_PassedParticleCombos(locSurvivedParticleCombos);

116

locParticleCombos.insert(locParticleCombos.end(), locSurvivedParticleCombos.begin(), locSurvivedParticleCombos.end());

117

}

118

119

dKinFitter.Reset_NewEvent();

120

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

121

{

122

const DParticleCombo* locParticleCombo = locParticleCombos[loc_i];

123

const DReaction* locReaction = locParticleCombo->Get_Reaction();

124

if(locReaction->Get_KinFitType() == d_NoFit)

125

continue; //don't do any kinematic fits!

126

127

if(dDebugLevel > 0)

128

cout << "Create kinematic fit constraints for event: " << eventnumber << endl;

129

130

map<const DKinFitParticle*, pair<Particle_t, deque<const DKinematicData*> > > locDecayingKinFitParticles;

131

deque<pair<DKinFitConstraint_VertexBase*, set<DKinFitConstraint_P4*> > > locSortedConstraints;

132

deque<DKinFitConstraint*> locOriginalConstraints;

133

if(!Create_KinFitConstraints(locParticleCombo, locDecayingKinFitParticles, locOriginalConstraints, locSortedConstraints))

134

continue; //sort-constraints failed: invalid! cannot setup kinfit

135

136

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

137

138

//check if previous combos would have resulted in duplicate results: if so, just copy the result (or abort fit if it failed last time)

139

//this could happen if you want to perform the exact same kinematic fit, but want to perform different DAnalysisActions on them

140

if(Handle_IfKinFitResultsWillBeIdentical(locParticleCombo, locOriginalConstraints, locEventRFBunch, locDecayingKinFitParticles))

141

continue; //identical results

142

143

//setup kinfit

144

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

145

if(!Setup_KinFit(locKinFitType, locOriginalConstraints, locEventRFBunch, locSortedConstraints))

146

{

147

if(dDebugLevel > 0)

148

cout << "Kinematic Fit Setup Failed." << endl;

149

deque<const DKinFitConstraint*> locConstOriginalConstraints(locOriginalConstraints.begin(), locOriginalConstraints.end());

150

dPreviouslyFailedFits.push_back(DPreviousFitInfo(locEventRFBunch, locConstOriginalConstraints, locDecayingKinFitParticles));

151

continue;

152

}

153

154

if(dDebugLevel > 0)

155

cout << "Perform Primary Kinematic Fit" << endl;

156

if(dKinFitter.Fit_Reaction()) //if fit fails: no kinfit results, BUT will still generate new DParticleCombo (using old info though!)

157

{

158

if(dDebugLevel > 0)

159

cout << "Kinematic Fit Converged for Event " << eventnumber << ", confidence = " << dKinFitter.Get_ConfidenceLevel() << endl;

160

Build_KinFitResults(locParticleCombo, locDecayingKinFitParticles, locOriginalConstraints);

161

}

162

else

163

{

164

deque<const DKinFitConstraint*> locConstOriginalConstraints(locOriginalConstraints.begin(), locOriginalConstraints.end());

165

dPreviouslyFailedFits.push_back(DPreviousFitInfo(locEventRFBunch, locConstOriginalConstraints, locDecayingKinFitParticles));

166

}

167

}

168

169

return NOERROR;

170

}

171

172

bool DKinFitResults_factory::Handle_IfKinFitResultsWillBeIdentical(const DParticleCombo* locParticleCombo, deque<DKinFitConstraint*> locConstraints_ToCheck, const DEventRFBunch* locRFBunch_ToCheck, map<const DKinFitParticle*, pair<Particle_t, deque<const DKinematicData*> > > locDecayingKinFitParticles_ToCheck)

173

{

174

#ifdef VTRACE

175

VT_TRACER("DKinFitResults_factory::Handle_IfKinFitResultsWillBeIdentical()");

176

#endif

177

//first check previously-passed results

178

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

179

{

180

deque<const DKinFitConstraint*> locPreviousOriginalConstraints;

181

_data[loc_j]->Get_OriginalKinFitConstraints(locPreviousOriginalConstraints);

182

183

set<const DParticleCombo*> locPreviousParticleCombos;

184

_data[loc_j]->Get_ParticleCombos(locPreviousParticleCombos);

185

186

//ok to just grab the rf bunch from the first combo: if they are different between combos, then they weren't needed anyway (else would be separate results)

187

const DEventRFBunch* locPreviousEventRFBunch = (*locPreviousParticleCombos.begin())->Get_EventRFBunch();

188

189

map<const DKinFitParticle*, pair<Particle_t, deque<const DKinematicData*> > > locDecayingKinFitParticles_CheckAgainst;

190

_data[loc_j]->Get_InputDecayingParticleInfo(locDecayingKinFitParticles_CheckAgainst);

191

192

if(!Check_IfKinFitResultsWillBeIdentical(locDecayingKinFitParticles_ToCheck, locDecayingKinFitParticles_CheckAgainst))

193

continue;

194

if(!Check_IfKinFitResultsWillBeIdentical(locConstraints_ToCheck, locPreviousOriginalConstraints, locRFBunch_ToCheck, locPreviousEventRFBunch))

195

continue;

196

if(dDebugLevel > 0)

197

cout << "kinfit results will be the same: register combo & skip fit" << endl;

198

199

//kinfit results will be the same: clone and save

200

_data[loc_j]->Add_ParticleCombo(locParticleCombo);

201

return true;

202

}

203

204

//now check previously-failed fits

205

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

206

{

207

if(!Check_IfKinFitResultsWillBeIdentical(locDecayingKinFitParticles_ToCheck, dPreviouslyFailedFits[loc_j].dDecayingParticles))

208

continue;

209

if(!Check_IfKinFitResultsWillBeIdentical(locConstraints_ToCheck, dPreviouslyFailedFits[loc_j].dOriginalConstraints, locRFBunch_ToCheck, dPreviouslyFailedFits[loc_j].dEventRFBunch))

210

continue;

211

if(dDebugLevel > 0)

212

cout << "kinfit results will be the same & will fail: skip fit" << endl;

213

return true;

214

}

215

return false;

216

}

217

218

bool DKinFitResults_factory::Check_IfKinFitResultsWillBeIdentical(map<const DKinFitParticle*, pair<Particle_t, deque<const DKinematicData*> > > locDecayingKinFitParticles_ToCheck, map<const DKinFitParticle*, pair<Particle_t, deque<const DKinematicData*> > > locDecayingKinFitParticles_CheckAgainst)

219

{

220

if(locDecayingKinFitParticles_ToCheck.size() != locDecayingKinFitParticles_CheckAgainst.size())

221

return false;

222

223

map<const DKinFitParticle*, pair<Particle_t, deque<const DKinematicData*> > >::iterator locIterator_ToCheck = locDecayingKinFitParticles_ToCheck.begin();

224

for(; locIterator_ToCheck != locDecayingKinFitParticles_ToCheck.end(); ++locIterator_ToCheck)

225

{

226

map<const DKinFitParticle*, pair<Particle_t, deque<const DKinematicData*> > >::iterator locIterator_CheckAgainst = locDecayingKinFitParticles_CheckAgainst.begin();

227

bool locMatchFoundFlag = false;

228

for(; locIterator_CheckAgainst != locDecayingKinFitParticles_CheckAgainst.end(); ++locIterator_CheckAgainst)

229

{

230

if(locIterator_ToCheck->second.first != locIterator_CheckAgainst->second.first)

231

continue;

232

deque<const DKinematicData*> locParticles_ToCheck = locIterator_ToCheck->second.second;

233

deque<const DKinematicData*> locParticles_CheckAgainst = locIterator_CheckAgainst->second.second;

234

if(locParticles_ToCheck.size() != locParticles_CheckAgainst.size())

235

continue;

236

237

//compare particles //note particles can be in a different order!

238

bool locAllParticlesMatchFlag = true;

239

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

240

{

241

bool locParticleMatchFoundFlag = false;

242

deque<const DKinematicData*>::iterator locDequeIterator = locParticles_CheckAgainst.begin();

243

for(; locDequeIterator != locParticles_CheckAgainst.end(); ++locDequeIterator)

244

{

245

if(locParticles_ToCheck[loc_i] != (*locDequeIterator))

246

continue;

247

locParticles_CheckAgainst.erase(locDequeIterator);

248

locParticleMatchFoundFlag = true;

249

break;

250

}

251

if(!locParticleMatchFoundFlag)

252

{

253

locAllParticlesMatchFlag = false;

254

break;

255

}

256

}

257

if(!locAllParticlesMatchFlag)

258

continue;

259

260

locMatchFoundFlag = true;

261

locDecayingKinFitParticles_CheckAgainst.erase(locIterator_CheckAgainst);

262

break;

263

}

264

if(!locMatchFoundFlag)

265

return false;

266

}

267

return true;

268

}

269

270

bool DKinFitResults_factory::Check_IfKinFitResultsWillBeIdentical(deque<DKinFitConstraint*> locConstraints_ToCheck, deque<const DKinFitConstraint*> locConstraints_CheckAgainst, const DEventRFBunch* locRFBunch_ToCheck, const DEventRFBunch* locRFBunch_CheckAgainst)

271

{

272

if(locConstraints_ToCheck.size() != locConstraints_CheckAgainst.size())

273

return false;

274

275

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

276

{

277

DKinFitConstraint_P4* locP4Constraint_ToCheck = dynamic_cast<DKinFitConstraint_P4*>(locConstraints_ToCheck[loc_i]);

278

DKinFitConstraint_Vertex* locVertexConstraint_ToCheck = dynamic_cast<DKinFitConstraint_Vertex*>(locConstraints_ToCheck[loc_i]);

279

DKinFitConstraint_Spacetime* locSpacetimeConstraint_ToCheck = dynamic_cast<DKinFitConstraint_Spacetime*>(locConstraints_ToCheck[loc_i]);

280

281

bool locMatchFoundFlag = false;

282

deque<const DKinFitConstraint*>::iterator locIterator = locConstraints_CheckAgainst.begin();

283

for(; locIterator != locConstraints_CheckAgainst.end(); ++locIterator)

284

{

285

if(locP4Constraint_ToCheck != NULL__null)

286

{

287

const DKinFitConstraint_P4* locP4Constraint_CheckAgainst = dynamic_cast<const DKinFitConstraint_P4*>(*locIterator);

288

if(locP4Constraint_CheckAgainst == NULL__null)

289

continue;

290

if(!Check_IfKinFitResultsWillBeIdentical(locP4Constraint_ToCheck, locP4Constraint_CheckAgainst))

291

continue;

292

locMatchFoundFlag = true;

293

locConstraints_CheckAgainst.erase(locIterator);

294

break;

295

}

296

if(locVertexConstraint_ToCheck != NULL__null)

297

{

298

const DKinFitConstraint_Vertex* locVertexConstraint_CheckAgainst = dynamic_cast<const DKinFitConstraint_Vertex*>(*locIterator);

299

if(locVertexConstraint_CheckAgainst == NULL__null)

300

continue;

301

if(!Check_IfKinFitResultsWillBeIdentical(static_cast<DKinFitConstraint_VertexBase*>(locVertexConstraint_ToCheck), static_cast<const DKinFitConstraint_VertexBase*>(locVertexConstraint_CheckAgainst)))

302

continue;

303

locMatchFoundFlag = true;

304

locConstraints_CheckAgainst.erase(locIterator);

305

break;

306

}

307

if(locSpacetimeConstraint_ToCheck != NULL__null)

308

{

309

const DKinFitConstraint_Spacetime* locSpacetimeConstraint_CheckAgainst = dynamic_cast<const DKinFitConstraint_Spacetime*>(*locIterator);

310

if(locSpacetimeConstraint_CheckAgainst == NULL__null)

311

continue;

312

if(!Check_IfKinFitResultsWillBeIdentical(static_cast<DKinFitConstraint_VertexBase*>(locSpacetimeConstraint_ToCheck), static_cast<const DKinFitConstraint_VertexBase*>(locSpacetimeConstraint_CheckAgainst)))

313

continue;

314

if(!Check_IfKinFitResultsWillBeIdentical(locSpacetimeConstraint_ToCheck, locSpacetimeConstraint_CheckAgainst, locRFBunch_ToCheck, locRFBunch_CheckAgainst))

315

continue;

316

locMatchFoundFlag = true;

317

locConstraints_CheckAgainst.erase(locIterator);

318

break;

319

}

320

}

321

if(!locMatchFoundFlag)

322

return false;

323

}

324

return true;

325

}

326

327

bool DKinFitResults_factory::Check_IfKinFitResultsWillBeIdentical(DKinFitConstraint_P4* locConstraint_ToCheck, const DKinFitConstraint_P4* locConstraint_CheckAgainst)

328

{

329

if(!Check_IfKinFitResultsWillBeIdentical(locConstraint_ToCheck->Get_InitialParticles(), locConstraint_CheckAgainst->Get_InitialParticles()))

330

return false;

331

if(!Check_IfKinFitResultsWillBeIdentical(locConstraint_ToCheck->Get_FinalParticles(), locConstraint_CheckAgainst->Get_FinalParticles()))

332

return false;

333

return true;

334

}

335

336

bool DKinFitResults_factory::Check_IfKinFitResultsWillBeIdentical(DKinFitConstraint_VertexBase* locConstraint_ToCheck, const DKinFitConstraint_VertexBase* locConstraint_CheckAgainst)

337

{

338

if(!Check_IfKinFitResultsWillBeIdentical(locConstraint_ToCheck->Get_FullConstrainParticles(), locConstraint_CheckAgainst->Get_FullConstrainParticles()))

339

return false;

340

if(!Check_IfKinFitResultsWillBeIdentical(locConstraint_ToCheck->Get_NoConstrainParticles(), locConstraint_CheckAgainst->Get_NoConstrainParticles()))

341

return false;

342

343

deque<pair<const DKinFitParticle*, bool> > locDecayingParticles_ToCheck = locConstraint_ToCheck->Get_DecayingParticles();

344

deque<pair<const DKinFitParticle*, bool> > locDecayingParticles_CheckAgainst = locConstraint_CheckAgainst->Get_DecayingParticles();

345

if(locDecayingParticles_ToCheck.size() != locDecayingParticles_CheckAgainst.size())

346

return false;

347

348

//compare sources of objects //note particles can be in a different order!

349

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

350

{

351

// get the input kinfit particle (particles in the constraints are clones!)

352

const DKinFitParticle* locDecayingParticle_ToCheck = dKinFitter.Get_InputKinFitParticle(locDecayingParticles_ToCheck[loc_i].first);

353

bool locMatchFoundFlag = false;

354

deque<pair<const DKinFitParticle*, bool> >::iterator locIterator = locDecayingParticles_CheckAgainst.begin();

355

for(; locIterator != locDecayingParticles_CheckAgainst.end(); ++locIterator)

356

{

357

// get the input kinfit particle (particles in the constraints are clones!)

358

const DKinFitParticle* locDecayingParticle_CheckAgainst = dKinFitter.Get_InputKinFitParticle((*locIterator).first);

359

360

//check pid & flag

361

if(locDecayingParticle_CheckAgainst->Get_PID() != locDecayingParticle_ToCheck->Get_PID())

362

continue;

363

if((*locIterator).second != locDecayingParticles_ToCheck[loc_i].second)

364

continue;

365

366

locDecayingParticles_CheckAgainst.erase(locIterator);

367

locMatchFoundFlag = true;

368

break;

369

}

370

if(!locMatchFoundFlag)

371

return false;

372

}

373

return true;

374

}

375

376

bool DKinFitResults_factory::Check_IfKinFitResultsWillBeIdentical(DKinFitConstraint_Spacetime* locConstraint_ToCheck, const DKinFitConstraint_Spacetime* locConstraint_CheckAgainst, const DEventRFBunch* locRFBunch_ToCheck, const DEventRFBunch* locRFBunch_CheckAgainst)

377

{

378

if(!Check_IfKinFitResultsWillBeIdentical(locConstraint_ToCheck->Get_OnlyConstrainTimeParticles(), locConstraint_CheckAgainst->Get_OnlyConstrainTimeParticles()))

379

return false;

380

if(locConstraint_ToCheck->Get_UseRFTimeFlag() != locConstraint_CheckAgainst->Get_UseRFTimeFlag())

381

return false;

382

if(locConstraint_ToCheck->Get_BeamParticle() != locConstraint_CheckAgainst->Get_BeamParticle())

383

return false;

384

if(locConstraint_ToCheck->Get_UseRFTimeFlag() && (locRFBunch_ToCheck != locRFBunch_CheckAgainst))

385

return false;

386

387

return true;

388

}

389

390

bool DKinFitResults_factory::Check_IfKinFitResultsWillBeIdentical(deque<const DKinFitParticle*> locParticles_ToCheck, deque<const DKinFitParticle*> locParticles_CheckAgainst)

391

{

392

if(locParticles_ToCheck.size() != locParticles_CheckAgainst.size())

Taking false branch

→

393

return false;

394

395

//compare sources of objects //note particles can be in a different order!

396

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

←

Loop condition is true. Entering loop body

→

397

{

398

// get the input kinfit particle (particles in the constraints are clones!)

399

const DKinFitParticle* locInputKinFitParticle_ToCheck = dKinFitter.Get_InputKinFitParticle(locParticles_ToCheck[loc_i]);

←

'locInputKinFitParticle_ToCheck' initialized here

→

400

bool locMatchFoundFlag = false;

401

deque<const DKinFitParticle*>::iterator locIterator = locParticles_CheckAgainst.begin();

402

for(; locIterator != locParticles_CheckAgainst.end(); ++locIterator)

←

Loop condition is true. Entering loop body

→

403

{

404

// get the input kinfit particle (particles in the constraints are clones!)

405

const DKinFitParticle* locInputKinFitParticle_CheckAgainst = dKinFitter.Get_InputKinFitParticle(*locIterator);

406

if(locInputKinFitParticle_ToCheck != locInputKinFitParticle_CheckAgainst)

←

Assuming 'locInputKinFitParticle_ToCheck' is equal to 'locInputKinFitParticle_CheckAgainst'

→

←

Taking false branch

→

407

continue;

408

if(locInputKinFitParticle_ToCheck == NULL__null)

←

Assuming pointer value is null

→

←

Taking true branch

→

409

{

410

//e.g. both particles missing/decaying/target: check type & pid

411

if(locInputKinFitParticle_CheckAgainst->Get_PID() != locInputKinFitParticle_ToCheck->Get_PID())

←

Called C++ object pointer is null

412

continue;

413

if(locInputKinFitParticle_CheckAgainst->Get_KinFitParticleType() != locInputKinFitParticle_ToCheck->Get_KinFitParticleType())

414

continue;

415

}

416

locParticles_CheckAgainst.erase(locIterator);

417

locMatchFoundFlag = true;

418

break;

419

}

420

if(!locMatchFoundFlag)

421

return false;

422

}

423

return true;

424

}

425

426

bool DKinFitResults_factory::Create_KinFitConstraints(const DParticleCombo* locParticleCombo, map<const DKinFitParticle*, pair<Particle_t, deque<const DKinematicData*> > >& locDecayingKinFitParticles, deque<DKinFitConstraint*>& locOriginalConstraints, deque<pair<DKinFitConstraint_VertexBase*, set<DKinFitConstraint_P4*> > >& locSortedConstraints)

427

{

428

#ifdef VTRACE

429

VT_TRACER("DKinFitResults_factory::Create_KinFitConstraints()");

430

#endif

431

if(dDebugLevel > 0)

432

cout << "DKinFitResults_factory: Create Particles" << endl;

433

434

locOriginalConstraints.clear();

435

locSortedConstraints.clear();

436

deque<deque<const DKinFitParticle*> > locInitialKinFitParticles, locFinalKinFitParticles;

437

438

const DParticleComboStep* locParticleComboStep;

439

int locDecayStepIndex;

440

const DReaction* locReaction = locParticleCombo->Get_Reaction();

441

const DKinematicData* locKinematicData;

442

const DKinFitParticle* locKinFitParticle;

443

Particle_t locPID;

444

DKinFitType locKinFitType = locReaction->Get_KinFitType();

445

446

//MAKE PARTICLES

447

locInitialKinFitParticles.clear();

448

locFinalKinFitParticles.clear();

449

map<size_t, const DKinFitParticle*> locDecayingParticleStepMap; //the map key is the step at which the particle decays at

450

bool locFirstParticleIsBeamFlag = (locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticleID() == Gamma);

451

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

452

{

453

deque<const DKinFitParticle*> locInitialKinFitParticles_Step;

454

deque<const DKinFitParticle*> locFinalKinFitParticles_Step;

455

locParticleComboStep = locParticleCombo->Get_ParticleComboStep(loc_i);

456

457

//initial particle

458

locPID = locParticleComboStep->Get_InitialParticleID();

459

if(locPID == Gamma)

460

{

461

const DBeamPhoton* locBeamPhoton = static_cast<const DBeamPhoton*>(locParticleComboStep->Get_InitialParticle());

462

locKinFitParticle = dKinFitter.Make_BeamParticle(locBeamPhoton);

463

locInitialKinFitParticles_Step.push_back(locKinFitParticle);

464

}

465

else if((!locFirstParticleIsBeamFlag) && (loc_i == 0)) //add decaying particle (production not specified, only decay)

466

{

467

locKinFitParticle = dKinFitter.Make_DecayingParticle(locPID);

468

deque<const DKinematicData*> locFinalDecayProducts;

469

locParticleCombo->Get_DecayChainParticles_Measured(loc_i, locFinalDecayProducts);

470

pair<Particle_t, deque<const DKinematicData*> > locDecayInfo(locPID, locFinalDecayProducts);

471

locDecayingKinFitParticles[locKinFitParticle] = locDecayInfo;

472

locInitialKinFitParticles_Step.push_back(locKinFitParticle);

473

locDecayingParticleStepMap[loc_i] = locKinFitParticle;

474

}

475

else //decaying particle, already added when it was a final particle

476

locInitialKinFitParticles_Step.push_back(locDecayingParticleStepMap[loc_i]);

477

478

//target particle

479

locPID = locParticleComboStep->Get_TargetParticleID();

480

if(locPID != Unknown)

481

{

482

locKinFitParticle = dKinFitter.Make_TargetParticle(locPID);

483

locInitialKinFitParticles_Step.push_back(locKinFitParticle);

484

}

485

486

//final state particles

487

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

488

{

489

locDecayStepIndex = locParticleComboStep->Get_DecayStepIndex(loc_j);

490

locKinematicData = locParticleComboStep->Get_FinalParticle(loc_j);

491

locPID = locParticleComboStep->Get_FinalParticleID(loc_j);

492

if(locDecayStepIndex == -1) //missing particle

493

{

494

locKinFitParticle = dKinFitter.Make_MissingParticle(locPID);

495

locFinalKinFitParticles_Step.push_back(locKinFitParticle);

496

}

497

else if(locDecayStepIndex >= 0) //add decaying particle

498

{

499

locKinFitParticle = dKinFitter.Make_DecayingParticle(locPID);

500

deque<const DKinematicData*> locFinalDecayProducts;

501

locParticleCombo->Get_DecayChainParticles_Measured(locDecayStepIndex, locFinalDecayProducts);

502

pair<Particle_t, deque<const DKinematicData*> > locDecayInfo(locPID, locFinalDecayProducts);

503

locDecayingKinFitParticles[locKinFitParticle] = locDecayInfo;

504

locFinalKinFitParticles_Step.push_back(locKinFitParticle);

505

locDecayingParticleStepMap[locDecayStepIndex] = locKinFitParticle;

506

}

507

else if(ParticleCharge(locKinematicData->PID()) == 0) //detected neutral

508

{

509

const DNeutralParticleHypothesis* locNeutralParticleHypothesis = static_cast<const DNeutralParticleHypothesis*>(locKinematicData);

510

if(locKinFitType == d_P4Fit) //make particle

511

locKinFitParticle = dKinFitter.Make_DetectedParticle(locNeutralParticleHypothesis);

512

else //make shower

513

locKinFitParticle = dKinFitter.Make_DetectedShower(locNeutralParticleHypothesis);

514

locFinalKinFitParticles_Step.push_back(locKinFitParticle);

515

}

516

else //detected charged track

517

{

518

const DChargedTrackHypothesis* locChargedTrackHypothesis = static_cast<const DChargedTrackHypothesis*>(locKinematicData);

519

locKinFitParticle = dKinFitter.Make_DetectedParticle(locChargedTrackHypothesis);

520

locFinalKinFitParticles_Step.push_back(locKinFitParticle);

521

}

522

}

523

524

//add particles to deques for constraint setup

525

locInitialKinFitParticles.push_back(locInitialKinFitParticles_Step);

526

locFinalKinFitParticles.push_back(locFinalKinFitParticles_Step);

527

}

528

if(dDebugLevel > 10)

529

cout << "DKinFitResults_factory: Particles Created." << endl;

530

531

//P4: Setup constraints

532

deque<DKinFitConstraint_P4*> locP4Constraints;

533

if((locKinFitType == d_P4Fit) || (locKinFitType == d_P4AndVertexFit) || (locKinFitType == d_P4AndSpacetimeFit))

534

{

535

if(dDebugLevel > 10)

536

cout << "DKinFitResults_factory: Setup P4 Constraints" << endl;

537

set<size_t> locStepIndicesToHandle;

538

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

539

locStepIndicesToHandle.insert(loc_i);

540

while(!locStepIndicesToHandle.empty())

541

{

542

deque<const DKinFitParticle*> locInitialKinFitParticles_P4, locFinalKinFitParticles_P4;

543

deque<size_t> locIncludedStepIndices;

544

545

bool locConstrainMassFlag = true; //note that the kinfitter ignores this flag when it doesn't apply (e.g. missing-mass/full-p4 constraint)

546

Setup_P4Constraint(locParticleCombo, (*locStepIndicesToHandle.begin()), locInitialKinFitParticles, locFinalKinFitParticles, locInitialKinFitParticles_P4, locFinalKinFitParticles_P4, locIncludedStepIndices, locConstrainMassFlag);

547

548

//remove steps included in the p4 constraint from the to-handle deque

549

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

550

locStepIndicesToHandle.erase(locIncludedStepIndices[loc_i]);

551

552

//make constraint

553

if(dDebugLevel > 10)

554

cout << "DKinFitResults_factory: Create P4 Constraint" << endl;

555

DKinFitConstraint_P4* locConstraint = dKinFitter.Make_P4Constraint(locInitialKinFitParticles_P4, locFinalKinFitParticles_P4, locConstrainMassFlag);

556

locOriginalConstraints.push_back(dynamic_cast<DKinFitConstraint*>(locConstraint));

557

locP4Constraints.push_back(locConstraint);

558

}

559

}

560

561

//VERTEX OR SPACETIME: Group particles by detached vertex (one deque for each constraint/vertex)

562

set<DKinFitConstraint_VertexBase*> locOriginalVertexBaseConstraints;

563

if((locKinFitType == d_VertexFit) || (locKinFitType == d_SpacetimeFit) || (locKinFitType == d_P4AndVertexFit) || (locKinFitType == d_P4AndSpacetimeFit))

564

{

565

if(dDebugLevel > 10)

566

cout << "DKinFitResults_factory: Setup Vertex & Spacetime Constraints" << endl;

567

set<size_t> locStepIndicesToHandle;

568

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

569

locStepIndicesToHandle.insert(loc_i);

570

while(!locStepIndicesToHandle.empty())

571

{

572

deque<const DKinFitParticle*> locInitialKinFitParticles_Vertex, locFinalKinFitParticles_Vertex;

573

deque<size_t> locIncludedStepIndices;

574

575

Setup_VertexConstraint(locParticleCombo, (*locStepIndicesToHandle.begin()), locInitialKinFitParticles, locFinalKinFitParticles, locInitialKinFitParticles_Vertex, locFinalKinFitParticles_Vertex, locIncludedStepIndices);

576

577

//remove steps included in the vertex constraint from the to-handle deque

578

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

579

locStepIndicesToHandle.erase(locIncludedStepIndices[loc_i]);

580

581

if((locKinFitType == d_VertexFit) || (locKinFitType == d_P4AndVertexFit))

582

{

583

DKinFitConstraint_Vertex* locConstraint = dKinFitter.Make_VertexConstraint(locInitialKinFitParticles_Vertex, locFinalKinFitParticles_Vertex, TVector3());

584

//string

585

string locConstraintString = string("#it{x}^{3}_{");

586

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

587

{

588

if((loc_j == 0) || (Is_FinalStateParticle(PDGtoPType(locInitialKinFitParticles_Vertex[loc_j]->Get_PID())) == 1))

589

locConstraintString += ParticleName_ROOT(PDGtoPType(locInitialKinFitParticles_Vertex[loc_j]->Get_PID()));

590

}

591

locConstraintString += "#rightarrow";

592

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

593

{

594

string locParticleString = ParticleName_ROOT(PDGtoPType(locFinalKinFitParticles_Vertex[loc_j]->Get_PID()));

595

if(locFinalKinFitParticles_Vertex[loc_j]->Get_KinFitParticleType() == d_MissingParticle)

596

locConstraintString += string("(") + locParticleString + string(")");

597

else

598

locConstraintString += locParticleString;

599

}

600

locConstraintString += "}";

601

locConstraint->Set_ConstraintString(locConstraintString);

602

//save

603

locOriginalConstraints.push_back(dynamic_cast<DKinFitConstraint*>(locConstraint));

604

locOriginalVertexBaseConstraints.insert(dynamic_cast<DKinFitConstraint_VertexBase*>(locConstraint));

605

}

606

if((locKinFitType == d_SpacetimeFit) || (locKinFitType == d_P4AndSpacetimeFit))

607

{

608

bool locUseRFTimeFlag = false;

609

if(locFirstParticleIsBeamFlag && (!locInitialKinFitParticles_Vertex.empty()))

610

locUseRFTimeFlag = ((locInitialKinFitParticles_Vertex[0]->Get_Charge() == 0) && (!(locInitialKinFitParticles_Vertex[0]->Get_Mass() > 0.0))); //true if photon

611

DKinFitConstraint_Spacetime* locConstraint = dKinFitter.Make_SpacetimeConstraint(locInitialKinFitParticles_Vertex, locFinalKinFitParticles_Vertex, locUseRFTimeFlag, TVector3(), 0.0);

612

//string

613

string locConstraintString = string("#it{x}^{4}_{");

614

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

615

{

616

if((loc_j == 0) || (Is_FinalStateParticle(PDGtoPType(locInitialKinFitParticles_Vertex[loc_j]->Get_PID())) == 1))

617

locConstraintString += ParticleName_ROOT(PDGtoPType(locInitialKinFitParticles_Vertex[loc_j]->Get_PID()));

618

}

619

locConstraintString += "#rightarrow";

620

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

621

{

622

string locParticleString = ParticleName_ROOT(PDGtoPType(locFinalKinFitParticles_Vertex[loc_j]->Get_PID()));

623

if(locFinalKinFitParticles_Vertex[loc_j]->Get_KinFitParticleType() == d_MissingParticle)

624

locConstraintString += string("(") + locParticleString + string(")");

625

else

626

locConstraintString += locParticleString;

627

}

628

locConstraintString += "}";

629

locConstraint->Set_ConstraintString(locConstraintString);

630

//save

631

locOriginalConstraints.push_back(dynamic_cast<DKinFitConstraint*>(locConstraint));

632

locOriginalVertexBaseConstraints.insert(dynamic_cast<DKinFitConstraint_VertexBase*>(locConstraint));

633

}

634

}

635

}

636

if(dDebugLevel > 10)

637

cout << "DKinFitResults_factory: All Constraints Created." << endl;

638

639

if(dDebugLevel > 10)

640

cout << "DKinFitResults_factory: Sort Constraints." << endl;

641

642

if(!dKinFitter.Sort_Constraints(locOriginalConstraints, locSortedConstraints))

643

return false; //invalid vertex constraints are skimmed, but invalid p4 constraints cause return-false!!

644

645

if(dDebugLevel > 10)

646

cout << "DKinFitResults_factory: Constraints Sorted." << endl;

647

648

//search for neutral showers whose vertex fits got rejected: replace them with neutral particles!

649

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

650

locOriginalVertexBaseConstraints.erase(locSortedConstraints[loc_i].first);

651

//only the rejected constraints remain in locOriginalVertexBaseConstraints. loop over them and look for neutrals

652

set<DKinFitConstraint_VertexBase*>::iterator locSetIterator = locOriginalVertexBaseConstraints.begin();

653

set<const DKinFitParticle*> locNeutralShowersToReplace;

654

for(; locSetIterator != locOriginalVertexBaseConstraints.end(); ++locSetIterator)

655

{

656

deque<const DKinFitParticle*> locNoConstrainParticles = (*locSetIterator)->Get_NoConstrainParticles();

657

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

658

{

659

if(locNoConstrainParticles[loc_i]->Get_IsNeutralShowerFlag()) //uh oh, need to create a particle for its p4 constraint

660

locNeutralShowersToReplace.insert(locNoConstrainParticles[loc_i]);

661

}

662

//check spacetime constraints

663

DKinFitConstraint_Spacetime* locSpacetimeConstraint = dynamic_cast<DKinFitConstraint_Spacetime*>(*locSetIterator);

664

if(locSpacetimeConstraint == NULL__null)

665

continue;

666

deque<const DKinFitParticle*> locOnlyConstrainTimeParticles = locSpacetimeConstraint->Get_OnlyConstrainTimeParticles();

667

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

668

{

669

if(locOnlyConstrainTimeParticles[loc_i]->Get_IsNeutralShowerFlag()) //uh oh, need to create a particle for its p4 constraint

670

locNeutralShowersToReplace.insert(locNoConstrainParticles[loc_i]);

671

}

672

}

673

if(!locNeutralShowersToReplace.empty())

674

{

675

//loop over p4 constraints, looking for the neutral showers we need to replace

676

map<const DKinFitParticle*, const DKinematicData*> locKinematicDataMapping;

677

dKinFitter.Get_ParticleMapping_InputToSource(locKinematicDataMapping);

678

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

679

{

680

deque<const DKinFitParticle*> locFinalParticles = locP4Constraints[loc_i]->Get_FinalParticles();

681

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

682

{

683

if(locNeutralShowersToReplace.find(locFinalParticles[loc_j]) == locNeutralShowersToReplace.end())

684

continue;

685

686

const DKinematicData* locKinematicData = locKinematicDataMapping[locFinalParticles[loc_j]];

687

const DNeutralParticleHypothesis* locNeutralParticleHypothesis = static_cast<const DNeutralParticleHypothesis*>(locKinematicData);

688

const DKinFitParticle* locNewKinFitParticle = dKinFitter.Make_DetectedParticle(locNeutralParticleHypothesis);

689

locP4Constraints[loc_i]->Replace_Particle(locFinalParticles[loc_j], false, locNewKinFitParticle); //replace!!!

690

}

691

}

692

}

693

694

//strip skimmed vertex/time constraints from original constraints

695

deque<DKinFitConstraint*>::iterator locConstraintIterator = locOriginalConstraints.begin();

696

for(; locConstraintIterator != locOriginalConstraints.end();)

697

{

698

if(dynamic_cast<DKinFitConstraint_VertexBase*>(*locConstraintIterator) == NULL__null)

699

{

700

++locConstraintIterator;

701

continue;

702

}

703

704

//see if can find this constraint in the sorted constraints (may have been skimmed out)

705

bool locConstraintFoundFlag = false;

706

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

707

{

708

if(locSortedConstraints[loc_j].first != (*locConstraintIterator))

709

continue;

710

locConstraintFoundFlag = true;

711

break;

712

}

713

if(locConstraintFoundFlag)

714

++locConstraintIterator;

715

else

716

locConstraintIterator = locOriginalConstraints.erase(locConstraintIterator); //skimmed out, don't include

717

}

718

719

return true;

720

}

721

722

bool DKinFitResults_factory::Setup_KinFit(DKinFitType locKinFitType, const deque<DKinFitConstraint*>& locOriginalConstraints, const DEventRFBunch* locEventRFBunch, deque<pair<DKinFitConstraint_VertexBase*, set<DKinFitConstraint_P4*> > >& locSortedConstraints)

723

{

724

#ifdef VTRACE

725

VT_TRACER("DKinFitResults_factory::Setup_KinFit()");

726

#endif

727

//group p4 constraints together

728

deque<DKinFitConstraint_P4*> locP4Constraints;

729

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

730

{

731

DKinFitConstraint_P4* locP4Constraint = dynamic_cast<DKinFitConstraint_P4*>(locOriginalConstraints[loc_i]);

732

if(locP4Constraint != NULL__null)

733

locP4Constraints.push_back(locP4Constraint);

734

}

735

736

//check if don't need to bother with initial vertex guesses: if so, just set constraints and return

737

if((locKinFitType == d_P4Fit) || locSortedConstraints.empty())

738

{

739

dKinFitter.Reset_NewFit();

740

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

741

dKinFitter.Set_Constraint(locP4Constraints[loc_i]);

742

return true;

743

}

744

745

double locRFTime = locEventRFBunch->dTime;

746

double locRFUncertainty = sqrt(locEventRFBunch->dTimeVariance);

747

748

//METHOD: for each vertex (in order):

749

//1) kinfit vertex only

750

//2) if this is the last vertex, then exit

751

//3) if possible (see 4), kinfit decay (vertex-p4) (all p4s at that vertex at once)

752

//this gives decaying particles a full covariance matrix so that they can be used in subsequent vertex kinfits

753

//4) if 3) is not possible because there are >= 2 locally-open-ended missing/decaying no-constrain particles (or missing particle with unknown pid): perform overall p4 fit (unless already done)

754

// (p4 only (mass constraints included where possible, all p4 constraints): get p4 cov matrices for all missing & decaying particles

755

bool locP4OnlyFitPerformedFlag = false;

756

map<const DKinFitParticle*, pair<const DKinFitParticle*, bool> > locReconstructedParticleMap; //map from orig particle to reconstructed particle //bool is fit good/bad (true/false)

757

758

//last resort: decaying particles that were reconstruced in a p4-only fit that were NOT at that current vertex

759

//if can't perform a vertex-p4 fit, use this info instead

760

map<const DKinFitParticle*, pair<const DKinFitParticle*, bool> > locReconstructedParticleLastResortMap; //last resort map from orig particle to reconstructed particle //bool is fit good/bad (true/false)

761

762

map<const DKinFitParticle*, const DKinFitParticle*> locReconstructedDetectedParticleMap; //map from orig decaying/missing particle to replacement "detected" particle

763

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

764

{

765

/************************************************************** VERTEX-ONLY FIT **************************************************************/

766

767

if(dDebugLevel > 10)

768

cout << "DKinFitResults_factory: Init Fit " << loc_i + 1 << " of " << locSortedConstraints.size() << ": Setup Vertex-only fit." << endl;

769

770

DKinFitConstraint_VertexBase* locOriginalVertexBaseConstraint = locSortedConstraints[loc_i].first;

771

DKinFitConstraint_Vertex* locOriginalVertexConstraint = dynamic_cast<DKinFitConstraint_Vertex*>(locOriginalVertexBaseConstraint);

772

DKinFitConstraint_Spacetime* locOriginalSpacetimeConstraint = dynamic_cast<DKinFitConstraint_Spacetime*>(locOriginalVertexBaseConstraint);

773

774

//may modify the original constraints for this fit

775

DKinFitConstraint_VertexBase* locThisFitVertexBaseConstraint = locOriginalVertexBaseConstraint;

776

DKinFitConstraint_Vertex* locThisFitVertexConstraint = locOriginalVertexConstraint;

777

DKinFitConstraint_Spacetime* locThisFitSpacetimeConstraint = locOriginalSpacetimeConstraint;

778

779

//treat previously-constrained decaying particles as detected in these constraints

780

bool locConstraintClonedFlag = false;

781

bool locPreviousFitFailedFlag = false;

782

deque<pair<const DKinFitParticle*, bool> > locDecayingParticles = locOriginalVertexBaseConstraint->Get_DecayingParticles();

783

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

784

{

785

map<const DKinFitParticle*, pair<const DKinFitParticle*, bool> >::iterator locMapIterator = locReconstructedParticleMap.find(locDecayingParticles[loc_j].first);

786

if(locMapIterator == locReconstructedParticleMap.end())

787

continue;

788

if(!locMapIterator->second.second)

789

{

790

//previous fit for this particle failed

791

if(locThisFitVertexBaseConstraint->Get_FullConstrainParticles().size() >= 2)

792

continue; //but this particle isn't needed for the vertex kinematic fit: don't contaminate

793

locPreviousFitFailedFlag = true; //particle is needed: don't kinfit, just return init vertex guess

794

//Note: the p4 of this particle is potentially from a failed p4 fit. It should be accurate enough for a crude guess though (and better than nothing)

795

}

796

if(dDebugLevel > 10)

797

cout << "DKinFitResults_factory: Reconstructed particle found, replace decaying particle." << endl;

798

//replace this decaying particle with a fully-constrained "detected" particle

799

if(!locConstraintClonedFlag) //don't touch the original! clone orig constraint without cloning particles

800

{

801

if(dDebugLevel > 15)

802

cout << "DKinFitResults_factory: Clone Fit." << endl;

803

if(locOriginalVertexConstraint != NULL__null)

804

{

805

locThisFitVertexConstraint = dKinFitter.Clone_KinFitConstraint_Vertex(locOriginalVertexConstraint);

806

locThisFitVertexBaseConstraint = dynamic_cast<DKinFitConstraint_VertexBase*>(locThisFitVertexConstraint);

807

}

808

else //spacetime

809

{

810

locThisFitSpacetimeConstraint = dKinFitter.Clone_KinFitConstraint_Spacetime(locOriginalSpacetimeConstraint);

811

locThisFitVertexBaseConstraint = dynamic_cast<DKinFitConstraint_VertexBase*>(locThisFitSpacetimeConstraint);

812

}

813

locConstraintClonedFlag = true;

814

}

815

//if not already "detected," create a new "detected" particle

816

const DKinFitParticle* locReconstructedParticle = locMapIterator->second.first;

817

const DKinFitParticle* locParticleToTreatAsDetected = locReconstructedParticle;

818

if(locReconstructedParticle->Get_KinFitParticleType() != d_DetectedParticle)

819

{

820

if(dDebugLevel > 15)

821

cout << "DKinFitResults_factory: Make new 'Detected' Particle." << endl;

822

locParticleToTreatAsDetected = dKinFitter.Make_DetectedParticle(locReconstructedParticle->Get_PID(), locReconstructedParticle->Get_Charge(), locReconstructedParticle->Get_Mass(), locReconstructedParticle->Get_SpacetimeVertex(), locReconstructedParticle->Get_Momentum(), locReconstructedParticle->Get_CovarianceMatrix());

823

}

824

locReconstructedDetectedParticleMap[locReconstructedParticle] = locParticleToTreatAsDetected;

825

//replace decaying particle with new "detected" one

826

locThisFitVertexBaseConstraint->Replace_DecayingParticle(locDecayingParticles[loc_j].first, locParticleToTreatAsDetected);

827

}

828

829

//do crude vertex guess: point on DOCA-line between the two closest (by doca) particles

830

deque<const DKinFitParticle*> locFullConstrainParticles = locThisFitVertexBaseConstraint->Get_FullConstrainParticles();

831

DVector3 locDVertexGuess = dAnalysisUtilities->Calc_CrudeVertex(locFullConstrainParticles);

832

TVector3 locVertexGuess(locDVertexGuess.X(), locDVertexGuess.Y(), locDVertexGuess.Z());

833

if(dDebugLevel > 20)

834

cout << "init vertex guess = " << locVertexGuess.X() << ", " << locVertexGuess.Y() << ", " << locVertexGuess.Z() << endl;

835

locThisFitVertexBaseConstraint->Set_VertexGuess(locVertexGuess);

836

locOriginalVertexBaseConstraint->Set_VertexGuess(locVertexGuess);

837

838

if((locFullConstrainParticles.size() < 2) || locPreviousFitFailedFlag)

839

continue; //you failed a previous setup fit earlier, and don't have the needed information. you'll have to live with what you've got

840

841

//setup fitter for vertex-only fit

842

dKinFitter.Reset_NewFit();

843

double locTimeGuess = 0.0;

844

const DKinFitConstraint_VertexBase* locInitFitVertexBaseConstraint = NULL__null;

845

if(locThisFitVertexConstraint != NULL__null)

846

{

847

dKinFitter.Set_Constraint(locThisFitVertexConstraint);

848

locInitFitVertexBaseConstraint = static_cast<const DKinFitConstraint_VertexBase*>(locThisFitVertexConstraint);

849

}

850

else

851

{

852

locTimeGuess = Calc_TimeGuess(locThisFitSpacetimeConstraint, locDVertexGuess, locRFTime);

853

locThisFitSpacetimeConstraint->Set_TimeGuess(locTimeGuess);

854

locOriginalSpacetimeConstraint->Set_TimeGuess(locTimeGuess);

855

dKinFitter.Set_Constraint(locThisFitSpacetimeConstraint);

856

locInitFitVertexBaseConstraint = static_cast<const DKinFitConstraint_VertexBase*>(locThisFitSpacetimeConstraint);

857

}

858

859

//perform fit

860

if(dDebugLevel > 0)

861

cout << "Perform init vertex guess kinematic fit" << endl;

862

bool locInitVertexFitStatus = dKinFitter.Fit_Reaction();

863

locInitFitVertexBaseConstraint = dynamic_cast<const DKinFitConstraint_VertexBase*>(dKinFitter.Get_OutputKinFitConstraint(locInitFitVertexBaseConstraint));

864

// it is not "necessary" for this fit to converge, because we already have a vertex/time guess (determined roughly above).

865

//thus don't return false here: if the fit fails (e.g. max #iterations), try to go on with the next iteration

866

if(locInitVertexFitStatus)

867

{

868

//save results: constraints were cloned during fit, so need to update init guesses for the next fit

869

locVertexGuess = locInitFitVertexBaseConstraint->Get_CommonVertex(); //else don't update guess!

870

if(dDebugLevel > 20)

871

cout << "init vertex guess = " << locVertexGuess.X() << ", " << locVertexGuess.Y() << ", " << locVertexGuess.Z() << endl;

872

locThisFitVertexBaseConstraint->Set_VertexGuess(locVertexGuess);

873

locOriginalVertexBaseConstraint->Set_VertexGuess(locVertexGuess);

874

if(locThisFitSpacetimeConstraint != NULL__null)

875

{

876

const DKinFitConstraint_Spacetime* locInitFitSpacetimeConstraint = dynamic_cast<const DKinFitConstraint_Spacetime*>(locInitFitVertexBaseConstraint);

877

locThisFitSpacetimeConstraint->Set_TimeGuess(locInitFitSpacetimeConstraint->Get_CommonTime());

878

locOriginalSpacetimeConstraint->Set_TimeGuess(locInitFitSpacetimeConstraint->Get_CommonTime());

879

}

880

}

881

882

/*************************************************************** CHECK P4 STATUS *************************************************************/

883

884

if(loc_i == (locSortedConstraints.size() - 1))

885

break; //have all vertex guesses, don't bother with p4

886

887

if((locKinFitType == d_VertexFit) || (locKinFitType == d_SpacetimeFit))

888

continue; //no p4 fits, and decaying particles in multiple vertex/time fits with no p4 fit is illegal (will be rejected later), so don't worry about it: continue

889

890

//see how many locally-open-ended no-constrain missing/decaying particles there are: start with decaying

891

//this is an init/intermediary fit; all previously-fully-defined decaying particles should be treated as detected, and should NOT be in the decaying vector

892

//note that fully-constrained decaying particles are treated as detected and thus aren't in the decaying particle vector

893

deque<const DKinFitParticle*> locLocallyOpenEndedNoConstrainParticles;

894

locDecayingParticles = locThisFitVertexConstraint->Get_DecayingParticles();

895

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

896

locLocallyOpenEndedNoConstrainParticles.push_back(locDecayingParticles[loc_i].first);

897

898

if(locLocallyOpenEndedNoConstrainParticles.empty())

899

continue; //nothing to gain by doing a vertex-p4 fit (no decaying particles to reconstruct that would be used in a future vertex constraint): just go to next fit

900

901

//see how many locally-open-ended no-constrain missing/decaying particles there are: get missing

902

deque<const DKinFitParticle*> locNoConstrainParticles = locThisFitVertexConstraint->Get_NoConstrainParticles();

903

const DKinFitParticle* locVertexConstrainedMissingParticle = NULL__null;

904

bool locUnknownMissingParticleFlag = false;

905

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

906

{

907

DKinFitParticleType locKinFitParticleType = locNoConstrainParticles[loc_j]->Get_KinFitParticleType();

908

if(locKinFitParticleType == d_MissingParticle)

909

{

910

if(locNoConstrainParticles[loc_j]->Get_PID() == 0)

911

locUnknownMissingParticleFlag = true;

912

locVertexConstrainedMissingParticle = dKinFitter.Get_OutputKinFitParticle(locNoConstrainParticles[loc_j]);

913

locLocallyOpenEndedNoConstrainParticles.push_back(locNoConstrainParticles[loc_j]);

914

}

915

}

916

917

/*************************************************************** VERTEX-P4 FIT ***************************************************************/

918

919

if((locLocallyOpenEndedNoConstrainParticles.size() < 2) && (!locUnknownMissingParticleFlag))

920

{

921

if(dDebugLevel > 10)

922

cout << "DKinFitResults_factory: Init Fit " << loc_i + 1 << " of " << locSortedConstraints.size() << ": Setup Vertex-P4 fit." << endl;

923

924

//perform a vertex-p4 fit to get: better vertex guess & full cov matrix of remaining missing/decaying particle

925

dKinFitter.Reset_NewFit();

926

const DKinFitConstraint_VertexBase* locP4VertexFitVertexBaseConstraint = NULL__null;

927

if(locThisFitVertexConstraint != NULL__null)

928

{

929

dKinFitter.Set_Constraint(locThisFitVertexConstraint);

930

locP4VertexFitVertexBaseConstraint = static_cast<const DKinFitConstraint_VertexBase*>(locThisFitVertexConstraint);

931

}

932

else

933

{

934

dKinFitter.Set_Constraint(locThisFitSpacetimeConstraint);

935

locP4VertexFitVertexBaseConstraint = static_cast<const DKinFitConstraint_VertexBase*>(locThisFitSpacetimeConstraint);

936

}

937

938

//get & potentially modify p4 constraints: treat previously-constrained decaying/missing particles as detected in these constraints

939

set<DKinFitConstraint_P4*> locTempP4Constraints = locSortedConstraints[loc_i].second;

940

for(set<DKinFitConstraint_P4*>::iterator locP4Iterator = locTempP4Constraints.begin(); locP4Iterator != locTempP4Constraints.end(); ++locP4Iterator)

941

{

942

DKinFitConstraint_P4* locOriginalP4Constraint = *locP4Iterator;

943

//get missing/decaying particles in this constraint

944

deque<pair<const DKinFitParticle*, bool> > locMissingOrDecayingParticles; //bool is true/false if in initial/final state

945

deque<const DKinFitParticle*> locInitialParticles = locOriginalP4Constraint->Get_InitialParticles();

946

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

947

{

948

DKinFitParticleType locKinFitParticleType = locInitialParticles[loc_j]->Get_KinFitParticleType();

949

if((locKinFitParticleType == d_DecayingParticle) || (locKinFitParticleType == d_MissingParticle))

950

locMissingOrDecayingParticles.push_back(pair<const DKinFitParticle*, bool>(locInitialParticles[loc_j], true));

951

}

952

deque<const DKinFitParticle*> locFinalParticles = locOriginalP4Constraint->Get_FinalParticles();

953

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

954

{

955

DKinFitParticleType locKinFitParticleType = locFinalParticles[loc_j]->Get_KinFitParticleType();

956

if((locKinFitParticleType == d_DecayingParticle) || (locKinFitParticleType == d_MissingParticle))

957

locMissingOrDecayingParticles.push_back(pair<const DKinFitParticle*, bool>(locFinalParticles[loc_j], false));

958

}

959

960

//loop over missing/decaying particles, replacing with reconstructed particles as available

961

DKinFitConstraint_P4* locClonedP4Constraint = NULL__null;

962

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

963

{

964

map<const DKinFitParticle*, pair<const DKinFitParticle*, bool> >::iterator locMapIterator = locReconstructedParticleMap.find(locMissingOrDecayingParticles[loc_j].first);

965

if(locMapIterator == locReconstructedParticleMap.end())

966

continue;

967

if(!locMapIterator->second.second)

968

{

969

locPreviousFitFailedFlag = true;

970

break; //this particle failed in a previous fit: cannot trust results. do not perform vertex-p4 fit

971

}

972

//replace this decaying/missing particle with a fully-constrained "detected" particle

973

974

//don't touch the original constraint! clone orig constraint (if haven't already) without cloning particles

975

if(locClonedP4Constraint == NULL__null)

976

locClonedP4Constraint = dKinFitter.Clone_KinFitConstraint_P4(locOriginalP4Constraint);

977

978

const DKinFitParticle* locParticleToTreatAsDetected = NULL__null;

979

map<const DKinFitParticle*, const DKinFitParticle*>::iterator locDetectedIterator = locReconstructedDetectedParticleMap.find(locMissingOrDecayingParticles[loc_j].first);

980

if(locDetectedIterator != locReconstructedDetectedParticleMap.end())

981

locParticleToTreatAsDetected = locDetectedIterator->second; //"detected" particle already made: use it

982

else

983

{

984

//create new "detected" particle

985

const DKinFitParticle* locReconstructedParticle = locMapIterator->second.first;

986

if(locInitVertexFitStatus)

987

{

988

if(locReconstructedParticle->Get_KinFitParticleType() == d_MissingParticle)

989

{

990

//combine info from above vertex fit & previous p3 fit of missing particle

991

//Note that if it is charged and in a B-field, then the missing momentum will be propagated between this input vertex position and the soon-to-be-kinfit one

992

//However, this distance is likely to be small, and the missing p4 will likely not change over this distance

993

TMatrixDSym locCovarianceMatrix = *(locVertexConstrainedMissingParticle->Get_CovarianceMatrix()); //contains v3 & t from current vertex fit

994

const TMatrixDSym* locP3CovarianceMatrix = locReconstructedParticle->Get_CovarianceMatrix(); //contains p3 from overall-p4 fit

995

//update with p3 from overall-p4 fit

996

for(unsigned int loc_q = 0; loc_q < 3; ++loc_q)

997

{

998

for(unsigned int loc_r = 0; loc_r < 3; ++loc_r)

999

locCovarianceMatrix(loc_q, loc_r) = (*locP3CovarianceMatrix)(loc_q, loc_r);

1000

}

1001

locParticleToTreatAsDetected = dKinFitter.Make_DetectedParticle(locReconstructedParticle->Get_PID(), locReconstructedParticle->Get_Charge(), locReconstructedParticle->Get_Mass(), locVertexConstrainedMissingParticle->Get_SpacetimeVertex(), locReconstructedParticle->Get_Momentum(), &locCovarianceMatrix);

1002

}

1003

else if(locReconstructedParticle->Get_KinFitParticleType() == d_DecayingParticle)

1004

1005

else //already reconstructed as detected, just copy it

1006

locParticleToTreatAsDetected = locReconstructedParticle;

1007

}

1008

else //cannot trust init vertex/time fit

1009

{

1010

TMatrixDSym locCovarianceMatrix = *(locReconstructedParticle->Get_CovarianceMatrix()); //contains p3 from overall-p4 fit

1011

//force ~infinite uncertainty on vertex params //it's only a setup fit anyway

1012

for(unsigned int loc_q = 0; loc_q < 4; ++loc_q)

1013

{

1014

for(unsigned int loc_r = 0; loc_r < 4; ++loc_r)

1015

locCovarianceMatrix(loc_q + 3, loc_r + 3) = 9.9E99;

1016

}

1017

1018

if(locReconstructedParticle->Get_KinFitParticleType() == d_MissingParticle)

1019

locParticleToTreatAsDetected = dKinFitter.Make_DetectedParticle(locReconstructedParticle->Get_PID(), locReconstructedParticle->Get_Charge(), locReconstructedParticle->Get_Mass(), TLorentzVector(locVertexGuess, locTimeGuess), locReconstructedParticle->Get_Momentum(), &locCovarianceMatrix);

1020

else if(locReconstructedParticle->Get_KinFitParticleType() == d_DecayingParticle)

1021

locParticleToTreatAsDetected = dKinFitter.Make_DetectedParticle(locReconstructedParticle->Get_PID(), locReconstructedParticle->Get_Charge(), locReconstructedParticle->Get_Mass(), TLorentzVector(locVertexGuess, locTimeGuess), locReconstructedParticle->Get_Momentum(), &locCovarianceMatrix);

1022

else //already reconstructed as detected, just copy it

1023

locParticleToTreatAsDetected = locReconstructedParticle;

1024

}

1025

locReconstructedDetectedParticleMap[locReconstructedParticle] = locParticleToTreatAsDetected;

1026

}

1027

1028

//replace missing/decaying particle with new "detected" one

1029

locClonedP4Constraint->Replace_Particle(locMissingOrDecayingParticles[loc_j].first, locMissingOrDecayingParticles[loc_j].second, locParticleToTreatAsDetected);

1030

}

1031

if(locPreviousFitFailedFlag)

1032

break;

1033

if(locClonedP4Constraint != NULL__null)

1034

dKinFitter.Set_Constraint(locClonedP4Constraint);

1035

else

1036

dKinFitter.Set_Constraint(locOriginalP4Constraint);

1037

}

1038

1039

//perform fit

1040

if(dDebugLevel > 0)

1041

cout << "Perform init p4-vertex guess kinematic fit" << endl;

1042

1043

bool locP4VertexFitStatus = (!locPreviousFitFailedFlag) ? dKinFitter.Fit_Reaction() : false;

1044

locP4VertexFitVertexBaseConstraint = dynamic_cast<const DKinFitConstraint_VertexBase*>(dKinFitter.Get_OutputKinFitConstraint(locP4VertexFitVertexBaseConstraint));

1045

if(locP4VertexFitStatus)

1046

{

1047

TVector3 locVertexGuess = locP4VertexFitVertexBaseConstraint->Get_CommonVertex();

1048

if(dDebugLevel > 20)

1049

cout << "init vertex guess = " << locVertexGuess.X() << ", " << locVertexGuess.Y() << ", " << locVertexGuess.Z() << endl;

1050

1051

//save v3/t results

1052

locOriginalVertexBaseConstraint->Set_VertexGuess(locVertexGuess);

1053

if(locOriginalSpacetimeConstraint != NULL__null)

1054

{

1055

const DKinFitConstraint_Spacetime* locInitFitSpacetimeConstraint = dynamic_cast<const DKinFitConstraint_Spacetime*>(locP4VertexFitVertexBaseConstraint);

1056

locOriginalSpacetimeConstraint->Set_TimeGuess(locInitFitSpacetimeConstraint->Get_CommonTime());

1057

}

1058

1059

//save reconstructed decaying/missing particles

1060

//assumes all particles reconstructed during vertex-p4 fit were in both p4 & vertex (as no-constrain particles) fits

1061

deque<const DKinFitParticle*> locNoConstrainParticles = locP4VertexFitVertexBaseConstraint->Get_NoConstrainParticles();

1062

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

1063

{

1064

DKinFitParticleType locKinFitParticleType = locNoConstrainParticles[loc_j]->Get_KinFitParticleType();

1065

if(locKinFitParticleType != d_DecayingParticle)

1066

continue;

1067

const DKinFitParticle* locOriginalKinFitParticle = dKinFitter.Get_InputKinFitParticle(locNoConstrainParticles[loc_j]);

1068

locReconstructedParticleMap[locOriginalKinFitParticle] = pair<const DKinFitParticle*, bool>(locNoConstrainParticles[loc_j], true);

1069

break;

1070

}

1071

continue;

1072

}

1073

//if vertex-p4 fit did not converge, move onto full-p4 fit to get decaying particle info

1074

}

1075

1076

/************************************************************* FULL, P4-ONLY FIT *************************************************************/

1077

1078

//too many unknown particles: cannot do local vertex-p4 fit OR vertex-p4 fit failed

1079

if(!locP4OnlyFitPerformedFlag)

1080

{

1081

//if not done already, perform a full, p4-only fit (mass constraints are applied as possible: especially important if inclusive fit)

1082

//save the resulting, reconstructed decaying/missing particles for use in later fits

1083

//will not use decaying particles reconstructed at other vertices, except as a last resort

1084

1085

if(dDebugLevel > 10)

1086

cout << "DKinFitResults_factory: Init Fit " << loc_i + 1 << " of " << locSortedConstraints.size() << ": Setup P4-Only fit." << endl;

1087

1088

//setup fit

1089

dKinFitter.Reset_NewFit();

1090

deque<const DKinFitConstraint_P4*> locInitFitP4Constraints;

1091

1092

//must treat any neutral showers as neutral particles (no vertex fit here)

1093

map<const DKinFitParticle*, const DKinematicData*> locKinematicDataMapping;

1094

dKinFitter.Get_ParticleMapping_InputToSource(locKinematicDataMapping);

1095

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

1096

{

1097

bool locConstraintClonedFlag = false;

1098

DKinFitConstraint_P4* locThisP4Constraint = locP4Constraints[loc_j];

1099

for(size_t loc_k = 0; loc_k < locP4Constraints[loc_j]->Get_FinalParticles().size(); ++loc_k)

1100

{

1101

const DKinFitParticle* locKinFitParticle = (locP4Constraints[loc_j]->Get_FinalParticles())[loc_k];

1102

if(!locKinFitParticle->Get_IsNeutralShowerFlag())

1103

continue;

1104

const DKinematicData* locKinematicData = locKinematicDataMapping[locKinFitParticle];

1105

const DNeutralParticleHypothesis* locNeutralParticleHypothesis = static_cast<const DNeutralParticleHypothesis*>(locKinematicData);

1106

const DKinFitParticle* locNewKinFitParticle = dKinFitter.Make_DetectedParticle(locNeutralParticleHypothesis);

1107

if(!locConstraintClonedFlag) //don't touch the original! clone orig constraint without cloning particles

1108

{

1109

if(dDebugLevel > 15)

1110

cout << "DKinFitResults_factory: Clone P4 Fit." << endl;

1111

locThisP4Constraint = dKinFitter.Clone_KinFitConstraint_P4(locThisP4Constraint);

1112

}

1113

locThisP4Constraint->Replace_Particle(locKinFitParticle, false, locNewKinFitParticle);

1114

}

1115

dKinFitter.Set_Constraint(locThisP4Constraint);

1116

locInitFitP4Constraints.push_back(locThisP4Constraint);

1117

}

1118

1119

//do fit

1120

if(dDebugLevel > 0)

1121

cout << "Perform init full-p4 kinematic fit" << endl;

1122

bool locP4OnlyFitStatus = dKinFitter.Fit_Reaction();

1123

1124

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

1125

locInitFitP4Constraints[loc_j] = dynamic_cast<const DKinFitConstraint_P4*>(dKinFitter.Get_OutputKinFitConstraint(locInitFitP4Constraints[loc_j]));

1126

1127

//save results: save all reconstructed decaying/missing particle info (if not saved already): p3 & cov p3

1128

//save the missing & decaying-no-constrain-particles-at-this-vertex in the reconstruction map

1129

//save decaying particles at other vertices TO A SPECIAL LAST RESORT MAP

1130

1131

//loop over created p4 constraints, grabbing the constrained particle (decaying/missing) in each one (and saving the results)

1132

deque<const DKinFitParticle*> locNoVertexConstrainParticles = locInitFitVertexBaseConstraint->Get_NoConstrainParticles();

1133

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

1134

{

1135

//even if the fit did not converge, locReconstructedParticle should have a semi-good guess for the p4

1136

const DKinFitParticle* locReconstructedParticle = locInitFitP4Constraints[loc_j]->Get_ConstrainedP4Particle();

1137

const DKinFitParticle* locOriginalParticle = dKinFitter.Get_InputKinFitParticle(locReconstructedParticle);

1138

if((locReconstructedParticle == NULL__null) || (locOriginalParticle == NULL__null))

1139

continue; //e.g. no missing particle

1140

//see if a reconstructed version of this particle has already been saved

1141

if(locReconstructedParticleMap.find(locOriginalParticle) != locReconstructedParticleMap.end())

1142

continue; //have already saved a reconstructed version of this particle, which is better than this new one. skip it

1143

1144

//if decaying particle defined at this vertex: save the v3 & v3 cov matrix

1145

const DKinFitParticle* locVertexFitResult = NULL__null;

1146

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

1147

{

1148

//need the original

1149

const DKinFitParticle* locVertexFitOriginal = dKinFitter.Get_InputKinFitParticle(locNoVertexConstrainParticles[loc_k]);

1150

if(locVertexFitOriginal == NULL__null)

1151

continue; //shouldn't be possible ...

1152

//compare originals

1153

if(locVertexFitOriginal != locOriginalParticle)

1154

continue;

1155

locVertexFitResult = locNoVertexConstrainParticles[loc_k];

1156

break;

1157

}

1158

1159

//save the result: missing particle

1160

if(locOriginalParticle->Get_KinFitParticleType() != d_DecayingParticle)

1161

{

1162

if(!locP4OnlyFitStatus)

1163

{

1164

//need a non-null covariance matrix

1165

TMatrixDSym locCovarianceMatrix(7, 7);

1166

locCovarianceMatrix.ResizeTo(7, 7);

1167

locReconstructedParticle = dKinFitter.Make_DetectedParticle(locReconstructedParticle->Get_PID(), locReconstructedParticle->Get_Charge(), locReconstructedParticle->Get_Mass(), locReconstructedParticle->Get_SpacetimeVertex(), locReconstructedParticle->Get_Momentum(), &locCovarianceMatrix);

1168

}

1169

locReconstructedParticleMap[locOriginalParticle] = pair<const DKinFitParticle*, bool>(locReconstructedParticle, locP4OnlyFitStatus);

1170

continue;

1171

}

1172

else if(locVertexFitResult == NULL__null) //decaying particle but not in this vertex fit

1173

{

1174

if(!locP4OnlyFitStatus)

1175

{

1176

//need a non-null covariance matrix

1177

TMatrixDSym locCovarianceMatrix(7, 7);

1178

locCovarianceMatrix.ResizeTo(7, 7);

1179

1180

}

1181

if(dDebugLevel > 10)

1182

cout << "DKinFitResults_factory: Save last resort particle for PID = " << locReconstructedParticle->Get_PID() << endl;

1183

locReconstructedParticleLastResortMap[locOriginalParticle] = pair<const DKinFitParticle*, bool>(locReconstructedParticle, locP4OnlyFitStatus);

1184

}

1185

else //decaying particle in this vertex fit

1186

{

1187

//combine info from above vertex fit with this p3 fit result: create new "detected" particle

1188

if(locInitVertexFitStatus)

1189

{

1190

//fit converged

1191

TMatrixDSym locCovarianceMatrix = *(locVertexFitResult->Get_CovarianceMatrix()); //contains v3 & t from current vertex fit

1192

if(locP4OnlyFitStatus)

1193

{

1194

const TMatrixDSym* locP3CovarianceMatrix = locReconstructedParticle->Get_CovarianceMatrix(); //contains p3 from overall-p4 fit

1195

//update with p3 from overall-p4 fit

1196

for(unsigned int loc_q = 0; loc_q < 3; ++loc_q)

1197

{

1198

for(unsigned int loc_r = 0; loc_r < 3; ++loc_r)

1199

locCovarianceMatrix(loc_q, loc_r) = (*locP3CovarianceMatrix)(loc_q, loc_r);

1200

}

1201

}

1202

const DKinFitParticle* locNewReconstructedParticle = dKinFitter.Make_DetectedParticle(locReconstructedParticle->Get_PID(), locReconstructedParticle->Get_Charge(), locReconstructedParticle->Get_Mass(), locVertexFitResult->Get_SpacetimeVertex(), locReconstructedParticle->Get_Momentum(), &locCovarianceMatrix);

1203

locReconstructedParticleMap[locOriginalParticle] = pair<const DKinFitParticle*, bool>(locNewReconstructedParticle, locP4OnlyFitStatus);

1204

}

1205

else

1206

{

1207

//cannot trust init vertex/time fit

1208

TMatrixDSym locCovarianceMatrix(7, 7);

1209

locCovarianceMatrix.ResizeTo(7, 7);

1210

//force ~infinite uncertainty on vertex params //it's only a setup fit anyway

1211

const DKinFitParticle* locNewReconstructedParticle = dKinFitter.Make_DetectedParticle(locReconstructedParticle->Get_PID(), locReconstructedParticle->Get_Charge(), locReconstructedParticle->Get_Mass(), TLorentzVector(locVertexGuess, locTimeGuess), locReconstructedParticle->Get_Momentum(), &locCovarianceMatrix);

1212

locReconstructedParticleMap[locOriginalParticle] = pair<const DKinFitParticle*, bool>(locNewReconstructedParticle, false); //fit failed

1213

}

1214

}

1215

}

1216

locP4OnlyFitPerformedFlag = true;

1217

continue;

1218

}

1219

1220

/**************************************************** USE LAST-RESORT P4-ONLY-FIT RESULTS ****************************************************/

1221

1222

//p4-only fit already performed, so instead, grab info for previously-fit decaying particles from last-resort reconstruction map

1223

//they weren't used earlier, because we were hoping to do a vertex-p4 fit to get those results

1224

//get needed decaying particles from last resort map, update their info with the latest vertex, push to the reconstruction map, and continue (don't do any fits)

1225

1226

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

1227

{

1228

//get from last-resort map

1229

const DKinFitParticle* locOriginalParticle = locLocallyOpenEndedNoConstrainParticles[loc_j];

1230

if(dDebugLevel > 10)

1231

cout << "DKinFitResults_factory: Get last resort particle for PID = " << locOriginalParticle->Get_PID() << endl;

1232

1233

map<const DKinFitParticle*, pair<const DKinFitParticle*, bool> >::iterator locLastResortIterator = locReconstructedParticleLastResortMap.find(locOriginalParticle);

1234

if(locLastResortIterator == locReconstructedParticleLastResortMap.end())

1235

continue; //e.g. missing

1236

1237

//combine info from above vertex fit with this p3 fit result: create new "detected" particle

1238

const DKinFitParticle* locLastResortParticle = locLastResortIterator->second.first;

1239

if(!locInitVertexFitStatus)

1240

{

1241

//can't improve, use as is

1242

locReconstructedParticleMap[locOriginalParticle] = pair<const DKinFitParticle*, bool>(locLastResortParticle, false);

1243

continue;

1244

}

1245

1246

const DKinFitParticle* locVertexFitResult = dKinFitter.Get_OutputKinFitParticle(locOriginalParticle);

1247

if(!locLastResortIterator->second.second)

1248

{

1249

//p4 fit failed, but can at least set vertex position

1250

const DKinFitParticle* locNewReconstructedParticle = dKinFitter.Make_DetectedParticle(locLastResortParticle->Get_PID(), locLastResortParticle->Get_Charge(), locLastResortParticle->Get_Mass(), locVertexFitResult->Get_SpacetimeVertex(), locLastResortParticle->Get_Momentum(), locVertexFitResult->Get_CovarianceMatrix());

1251

locReconstructedParticleMap[locOriginalParticle] = pair<const DKinFitParticle*, bool>(locNewReconstructedParticle, false);

1252

continue;

1253

}

1254

//p4 & vertex fit successful: combine info into new particle

1255

TMatrixDSym locCovarianceMatrix = *(locVertexFitResult->Get_CovarianceMatrix()); //contains v3 & t from current vertex fit

1256

const TMatrixDSym* locP3CovarianceMatrix = locLastResortParticle->Get_CovarianceMatrix(); //contains p3 from overall-p4 fit

1257

//update with p3 from overall-p4 fit

1258

for(unsigned int loc_q = 0; loc_q < 3; ++loc_q)

1259

{

1260

for(unsigned int loc_r = 0; loc_r < 3; ++loc_r)

1261

locCovarianceMatrix(loc_q, loc_r) = (*locP3CovarianceMatrix)(loc_q, loc_r);

1262

}

1263

1264

locReconstructedParticleMap[locOriginalParticle] = pair<const DKinFitParticle*, bool>(locNewReconstructedParticle, true);

1265

}

1266

}

1267

1268

if(dDebugLevel > 10)

1269

cout << "DKinFitResults_factory: Init Guesses found; set all constraints for master fit." << endl;

1270

1271

dKinFitter.Reset_NewFit();

1272

1273

//ADD CONSTRAINTS

1274

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

1275

dKinFitter.Set_Constraint(locP4Constraints[loc_i]);

1276

1277

//use locOriginalConstraints to keep constraints in same order (doesn't really matter, but confidence level histogram looks better this way)

1278

const DKinFitParticle* locBeamParticleForRF = NULL__null;

1279

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

1280

{

1281

if(dynamic_cast<DKinFitConstraint_VertexBase*>(locOriginalConstraints[loc_i]) == NULL__null)

1282

continue;

1283

1284

DKinFitConstraint_Vertex* locVertexConstraint = dynamic_cast<DKinFitConstraint_Vertex*>(locOriginalConstraints[loc_i]);

1285

if(locVertexConstraint != NULL__null)

1286

{

1287

dKinFitter.Set_Constraint(locVertexConstraint);

1288

if(dDebugLevel > 10)

1289

{

1290

TVector3 locVertex = locVertexConstraint->Get_CommonVertex();

1291

cout << "final vertex guess (xyz) = " << locVertex.X() << ", " << locVertex.Y() << ", " << locVertex.Z() << endl;

1292

}

1293

continue;

1294

}

1295

DKinFitConstraint_Spacetime* locSpacetimeConstraint = dynamic_cast<DKinFitConstraint_Spacetime*>(locOriginalConstraints[loc_i]);

1296

if(locSpacetimeConstraint != NULL__null)

1297

{

1298

if((locBeamParticleForRF == NULL__null) && locSpacetimeConstraint->Get_UseRFTimeFlag())

1299

locBeamParticleForRF = locSpacetimeConstraint->Get_BeamParticle();

1300

dKinFitter.Set_Constraint(locSpacetimeConstraint);

1301

if(dDebugLevel > 10)

1302

{

1303

TLorentzVector locLorentzVector = locSpacetimeConstraint->Get_CommonSpacetimeVertex();

1304

cout << "final spacetime guess (xyzt) = " << locLorentzVector.X() << ", " << locLorentzVector.Y() << ", " << locLorentzVector.Z() << ", " << locLorentzVector.T() << endl;

1305

}

1306

}

1307

}

1308

1309

//Add RF time if needed

1310

if(locBeamParticleForRF != NULL__null)

1311

dKinFitter.Set_RFTime(locRFTime, locRFUncertainty, locBeamParticleForRF);

1312

1313

return true;

1314

}

1315

1316

void DKinFitResults_factory::Setup_P4Constraint(const DParticleCombo* locParticleCombo, size_t locStepIndex, const deque<deque<const DKinFitParticle*> >& locInitialKinFitParticles, const deque<deque<const DKinFitParticle*> >& locFinalKinFitParticles, deque<const DKinFitParticle*>& locInitialKinFitParticles_P4, deque<const DKinFitParticle*>& locFinalKinFitParticles_P4, deque<size_t>& locIncludedStepIndices, bool& locConstrainMassFlag)

1317

{

1318

locIncludedStepIndices.push_back(locStepIndex);

1319

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

1320

Particle_t locPID;

1321

1322

//initial particle

1323

locPID = locParticleComboStep->Get_InitialParticleID();

1324

if(locPID == Gamma)

1325

locInitialKinFitParticles_P4.push_back(locInitialKinFitParticles[locStepIndex][0]);

1326

else //decaying particle

1327

{

1328

locConstrainMassFlag = locParticleCombo->Get_ApplyKinFitMassConstraintOnInitialParticleFlag(locStepIndex);

1329

if(IsFixedMass(locPID)) //else is decaying resonance particle: don't constrain (don't recurse either (already done!))

1330

locInitialKinFitParticles_P4.push_back(locInitialKinFitParticles[locStepIndex][0]);

1331

}

1332

1333

//target particle

1334

if(locParticleComboStep->Get_TargetParticleID() != Unknown)

1335

locInitialKinFitParticles_P4.push_back(locInitialKinFitParticles[locStepIndex][1]);

1336

1337

//final state particles

1338

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

1339

{

1340

int locDecayStepIndex = locParticleComboStep->Get_DecayStepIndex(loc_j);

1341

locPID = locParticleComboStep->Get_FinalParticleID(loc_j);

1342

if(locDecayStepIndex == -1) //missing particle

1343

locFinalKinFitParticles_P4.push_back(locFinalKinFitParticles[locStepIndex][loc_j]);

1344

else if(locDecayStepIndex >= 0) //decaying particle

1345

{

1346

if(!IsFixedMass(locPID)) //go to the next step!!

1347

Setup_P4Constraint(locParticleCombo, locDecayStepIndex, locInitialKinFitParticles, locFinalKinFitParticles, locInitialKinFitParticles_P4, locFinalKinFitParticles_P4, locIncludedStepIndices, locConstrainMassFlag);

1348

else

1349

locFinalKinFitParticles_P4.push_back(locFinalKinFitParticles[locStepIndex][loc_j]);

1350

}

1351

else //detected particle or shower

1352

locFinalKinFitParticles_P4.push_back(locFinalKinFitParticles[locStepIndex][loc_j]);

1353

}

1354

}

1355

1356

void DKinFitResults_factory::Setup_VertexConstraint(const DParticleCombo* locParticleCombo, size_t locStepIndex, const deque<deque<const DKinFitParticle*> >& locInitialKinFitParticles, const deque<deque<const DKinFitParticle*> >& locFinalKinFitParticles, deque<const DKinFitParticle*>& locInitialKinFitParticles_Vertex, deque<const DKinFitParticle*>& locFinalKinFitParticles_Vertex, deque<size_t>& locIncludedStepIndices)

1357

{

1358

locIncludedStepIndices.push_back(locStepIndex);

1359

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

1360

Particle_t locPID;

1361

1362

//initial particle

1363

locPID = locParticleComboStep->Get_InitialParticleID();

1364

if(locPID == Gamma)

1365

locInitialKinFitParticles_Vertex.push_back(locInitialKinFitParticles[locStepIndex][0]);

1366

else //decaying particle

1367

locInitialKinFitParticles_Vertex.push_back(locInitialKinFitParticles[locStepIndex][0]);

1368

1369

//target particle

1370

if(locParticleComboStep->Get_TargetParticleID() != Unknown)

1371

locInitialKinFitParticles_Vertex.push_back(locInitialKinFitParticles[locStepIndex][1]);

1372

1373

//final state particles

1374

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

1375

{

1376

int locDecayStepIndex = locParticleComboStep->Get_DecayStepIndex(loc_j);

1377

locPID = locParticleComboStep->Get_FinalParticleID(loc_j);

1378

if(locDecayStepIndex == -1) //missing particle

1379

locFinalKinFitParticles_Vertex.push_back(locFinalKinFitParticles[locStepIndex][loc_j]);

1380

else if(locDecayStepIndex >= 0) //decaying particle

1381

{

1382

if(IsDetachedVertex(locPID))

1383

{

1384

if(dLinkVerticesFlag)

1385

locFinalKinFitParticles_Vertex.push_back(locFinalKinFitParticles[locStepIndex][loc_j]);

1386

}

1387

else //go to the next step!!

1388

Setup_VertexConstraint(locParticleCombo, locDecayStepIndex, locInitialKinFitParticles, locFinalKinFitParticles, locInitialKinFitParticles_Vertex, locFinalKinFitParticles_Vertex, locIncludedStepIndices);

1389

}

1390

else //detected particle or shower

1391

locFinalKinFitParticles_Vertex.push_back(locFinalKinFitParticles[locStepIndex][loc_j]);

1392

}

1393

}

1394

1395

double DKinFitResults_factory::Calc_TimeGuess(const DKinFitConstraint_Spacetime* locConstraint, DVector3 locVertexGuess, double locRFTime)

1396

{

1397

//if RF: propagate RF time to the vertex-z of the init vertex guess

1398

bool locUseRFTimeFlag = locConstraint->Get_UseRFTimeFlag();

1399

if(locUseRFTimeFlag)

1400

return (locRFTime + (locVertexGuess.Z() - dTargetZCenter)/29.9792458);

1401

1402

//propagate each track time to the DOCA to the init vertex guess and average them

1403

deque<const DKinFitParticle*> locTimeFindParticles = locConstraint->Get_FullConstrainParticles();

1404

deque<const DKinFitParticle*> locOnlyTimeFindParticles = locConstraint->Get_OnlyConstrainTimeParticles();

1405

locTimeFindParticles.insert(locTimeFindParticles.end(), locOnlyTimeFindParticles.begin(), locOnlyTimeFindParticles.end());

1406

return dAnalysisUtilities->Calc_CrudeTime(locTimeFindParticles, locVertexGuess);

1407

}

1408

1409

void DKinFitResults_factory::Build_KinFitResults(const DParticleCombo* locParticleCombo, const map<const DKinFitParticle*, pair<Particle_t, deque<const DKinematicData*> > >& locInitDecayingKinFitParticles, deque<DKinFitConstraint*>& locOriginalConstraints)

1410

{

1411

const DReaction* locReaction = locParticleCombo->Get_Reaction();

1412

DKinFitResults* locKinFitResults = new DKinFitResults();

1413

locKinFitResults->Add_ParticleCombo(locParticleCombo);

1414

1415

locKinFitResults->Set_KinFitType(locReaction->Get_KinFitType());

1416

locKinFitResults->Set_ConfidenceLevel(dKinFitter.Get_ConfidenceLevel());

1417

locKinFitResults->Set_ChiSq(dKinFitter.Get_ChiSq());

1418

locKinFitResults->Set_NDF(dKinFitter.Get_NDF());

1419

locKinFitResults->Set_VEta(dKinFitter.Get_VEta());

1420

locKinFitResults->Set_VXi(dKinFitter.Get_VXi());

1421

locKinFitResults->Set_V(dKinFitter.Get_V());

1422

locKinFitResults->Set_NumConstraints(dKinFitter.Get_NumConstraintEquations());

1423

locKinFitResults->Set_NumUnknowns(dKinFitter.Get_NumUnknowns());

1424

locKinFitResults->Set_OriginalKinFitConstraints(locOriginalConstraints);

1425

locKinFitResults->Set_InputDecayingParticleInfo(locInitDecayingKinFitParticles);

1426

1427

//save the constraints

1428

deque<const DKinFitConstraint*> locKinFitConstraints;

1429

dKinFitter.Get_KinFitConstraints(locKinFitConstraints);

1430

1431

string locBeamTargetString = "";

1432

if(locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticleID() == Gamma)

1433

locBeamTargetString += ParticleName_ROOT(Gamma);

1434

if(locParticleCombo->Get_ParticleComboStep(0)->Get_TargetParticleID() != Unknown)

1435

locBeamTargetString += ParticleName_ROOT(locParticleCombo->Get_ParticleComboStep(0)->Get_TargetParticleID());

1436

locBeamTargetString += "#rightarrow";

1437

1438

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

1439

{

1440

const DKinFitConstraint_P4* locP4Constraint = dynamic_cast<const DKinFitConstraint_P4*>(locKinFitConstraints[loc_i]);

1441

if(locP4Constraint != NULL__null)

1442

{

1443

string locConstraintString = "";

1444

if(locP4Constraint->Get_IsActualP4ConstraintFlag())

1445

locConstraintString = "#it{p}^{4}";

1446

else if(locP4Constraint->Get_ConstrainInitialParticleMassFlag())

1447

{

1448

const DKinFitParticle* locKinFitParticle = (locP4Constraint->Get_InitialParticles())[0];

1449

locConstraintString = string("#it{m}_{") + string(ParticleName_ROOT(PDGtoPType(locKinFitParticle->Get_PID()))) + string("}");

1450

}

1451

(const_cast<DKinFitConstraint_P4*>(locP4Constraint))->Set_ConstraintString(locConstraintString);

1452

}

1453

else //vertex/time constraint, set colors based on if used to constrain

1454

{

1455

const DKinFitConstraint_VertexBase* locVertexBaseConstraint = dynamic_cast<const DKinFitConstraint_VertexBase*>(locKinFitConstraints[loc_i]);

1456

TString locConstraintTString = locVertexBaseConstraint->Get_ConstraintString();

1457

deque<const DKinFitParticle*> locFullConstrainParticles = locVertexBaseConstraint->Get_FullConstrainParticles();

1458

set<int> locPIDsReplaced;

1459

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

1460

{

1461

int locPID = locFullConstrainParticles[loc_j]->Get_PID();

1462

if(locFullConstrainParticles[loc_j]->Get_KinFitParticleType() == d_BeamParticle)

1463

{

1464

string locParticleString = ParticleName_ROOT(PDGtoPType(locPID));

1465

locConstraintTString.Replace(0, locParticleString.size(), (string("#color[4]{") + locParticleString + string("}")).c_str());

1466

locBeamTargetString.insert(locParticleString.size(), "}");

1467

continue;

1468

}

1469

if(locPIDsReplaced.find(locPID) != locPIDsReplaced.end())

1470

continue;

1471

string locParticleString = ParticleName_ROOT(PDGtoPType(locPID));

1472

string locMissingParticleString = string("(") + locParticleString + string(")");

1473

locConstraintTString.ReplaceAll(locMissingParticleString.c_str(), "__TEMPORARY__"); //temporary //don't want to make missing particles blue!

1474

locConstraintTString.ReplaceAll(locBeamTargetString.c_str(), "__BEAMTARGET__"); //temporary //don't want to make beam/target particles blue! (here)

1475

1476

string locNewParticleString = string("#color[4]{") + locParticleString + string("}"); //blue

1477

locConstraintTString.ReplaceAll(locParticleString.c_str(), locNewParticleString.c_str());

1478

locConstraintTString.ReplaceAll("__TEMPORARY__", locMissingParticleString.c_str()); //restore

1479

locConstraintTString.ReplaceAll("__BEAMTARGET__", locBeamTargetString.c_str()); //restore

1480

locPIDsReplaced.insert(locPID);

1481

}

1482

const DKinFitConstraint_Spacetime* locSpacetimeConstraint = dynamic_cast<const DKinFitConstraint_Spacetime*>(locVertexBaseConstraint);

1483

if(locSpacetimeConstraint != NULL__null)

1484

{

1485

deque<const DKinFitParticle*> locOnlyConstrainTimeParticles = locSpacetimeConstraint->Get_OnlyConstrainTimeParticles();

1486

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

1487

{

1488

int locPID = locOnlyConstrainTimeParticles[loc_j]->Get_PID();

1489

if(locOnlyConstrainTimeParticles[loc_j]->Get_KinFitParticleType() == d_BeamParticle)

1490

{

1491

string locParticleString = ParticleName_ROOT(PDGtoPType(locPID));

1492

locConstraintTString.Replace(0, locParticleString.size(), (string("#color[4]{") + locParticleString + string("}")).c_str());

1493

locBeamTargetString.insert(locParticleString.size(), "}");

1494

continue;

1495

}

1496

if(locPIDsReplaced.find(locPID) != locPIDsReplaced.end())

1497

continue;

1498

string locParticleString = ParticleName_ROOT(PDGtoPType(locPID));

1499

string locMissingParticleString = string("(") + locParticleString + string(")");

1500

locConstraintTString.ReplaceAll(locMissingParticleString.c_str(), "__TEMPORARY__"); //temporary //don't want to make missing particles blue!

1501

locConstraintTString.ReplaceAll(locBeamTargetString.c_str(), "__BEAMTARGET__"); //temporary //don't want to make beam/target particles blue! (here)

1502

1503

string locNewParticleString = string("#color[4]{") + locParticleString + string("}"); //blue

1504

locConstraintTString.ReplaceAll(locParticleString.c_str(), locNewParticleString.c_str());

1505

locConstraintTString.ReplaceAll("__TEMPORARY__", locMissingParticleString.c_str()); //restore

1506

locConstraintTString.ReplaceAll("__BEAMTARGET__", locBeamTargetString.c_str()); //restore

1507

locPIDsReplaced.insert(locPID);

1508

}

1509

}

1510

(const_cast<DKinFitConstraint_VertexBase*>(locVertexBaseConstraint))->Set_ConstraintString((const char*)locConstraintTString);

1511

}

1512

}

1513

1514

locKinFitResults->Set_KinFitConstraints(locKinFitConstraints);

1515

1516

//map of particle data

1517

map<const DKinFitParticle*, const DKinematicData*> locParticleMapping_Output;

1518

dKinFitter.Get_ParticleMapping_OutputToSource(locParticleMapping_Output);

1519

locKinFitResults->Set_ParticleMapping(locParticleMapping_Output);

1520

1521

//reverse map of particle data (& missing)

1522

map<const DKinematicData*, const DKinFitParticle*> locReverseParticleMapping;

1523

for(map<const DKinFitParticle*, const DKinematicData*>::iterator locIterator = locParticleMapping_Output.begin(); locIterator != locParticleMapping_Output.end(); ++locIterator)

1524

{

1525

if(locIterator->first->Get_KinFitParticleType() == d_MissingParticle)

1526

locKinFitResults->Set_MissingParticle(locIterator->first); //missing

1527

if(locIterator->second == NULL__null)

1528

continue; //not detected (e.g. missing, decaying, target)

1529

locReverseParticleMapping[locIterator->second] = locIterator->first;

1530

}

1531

locKinFitResults->Set_ReverseParticleMapping(locReverseParticleMapping);

1532

1533

//decaying particles //input to function is the constructed decaying particles; must save the final decaying particles

1534

deque<const DKinFitParticle*> locKinFitParticles;

1535

dKinFitter.Get_KinFitParticles(locKinFitParticles);

1536

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

1537

{

1538

if(locKinFitParticles[loc_i]->Get_KinFitParticleType() != d_DecayingParticle)

1539

continue;

1540

const DKinFitParticle* locInputKinFitParticle = dKinFitter.Get_InputKinFitParticle(locKinFitParticles[loc_i]);

1541

map<const DKinFitParticle*, pair<Particle_t, deque<const DKinematicData*> > >::const_iterator locDecayingIterator;

1542

locDecayingIterator = locInitDecayingKinFitParticles.find(locInputKinFitParticle);

1543

locKinFitResults->Add_DecayingParticle(locDecayingIterator->second, locKinFitParticles[loc_i]);

1544

}

1545

1546

//pulls

1547

map<const DKinematicData*, map<DKinFitPullType, double> > locPulls;

1548

dKinFitter.Get_Pulls(locPulls);

1549

locKinFitResults->Set_Pulls(locPulls);

1550

1551

_data.push_back(locKinFitResults);

1552

}

1553

1554

//------------------

1555

// erun

1556

//------------------

1557

jerror_t DKinFitResults_factory::erun(void)

1558

{

1559

return NOERROR;

1560

}

1561

1562

//------------------

1563

// fini

1564

//------------------

1565

jerror_t DKinFitResults_factory::fini(void)

1566

{

1567

return NOERROR;

1568

}

1569