libraries/ANALYSIS/DEventWriterROOT.cc

Bug Summary

File:	libraries/ANALYSIS/DEventWriterROOT.cc
Location:	line 1177, column 70
Description:	Division by zero

Annotated Source Code

#include "DEventWriterROOT.h"

DEventWriterROOT::DEventWriterROOT(JEventLoop* locEventLoop)

{

dInitNumThrownArraySize = 20;

dInitNumBeamArraySize = 20;

dInitNumTrackArraySize = 50;

dInitNumShowerArraySize = 15;

dInitNumComboArraySize = 100;

//BEWARE: IF THIS IS CHANGED, CHANGE IN THE BLUEPRINT FACTORY AND THE ANALYSIS UTILITIES ALSO!!

dTrackSelectionTag = "PreSelect";

dShowerSelectionTag = "PreSelect";

gPARMS->SetDefaultParameter("COMBO:TRACK_SELECT_TAG", dTrackSelectionTag);

gPARMS->SetDefaultParameter("COMBO:SHOWER_SELECT_TAG", dShowerSelectionTag);

locEventLoop->GetSingle(dAnalysisUtilities);

vector<const DReaction*> locReactions;

Get_Reactions(locEventLoop, locReactions);

//CREATE & INITIALIZE ANALYSIS ACTIONS

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

{

if(!locReactions[loc_i]->Get_EnableTTreeOutputFlag())

continue;

dCutActionMap_ThrownTopology[locReactions[loc_i]] = new DCutAction_ThrownTopology(locReactions[loc_i], true);

dCutActionMap_ThrownTopology[locReactions[loc_i]]->Initialize(locEventLoop);

dCutActionMap_TrueCombo[locReactions[loc_i]] = new DCutAction_TrueCombo(locReactions[loc_i], 5.73303E-7, true); //+/- 5sigma

dCutActionMap_TrueCombo[locReactions[loc_i]]->Initialize(locEventLoop);

dCutActionMap_BDTSignalCombo[locReactions[loc_i]] = new DCutAction_BDTSignalCombo(locReactions[loc_i], 5.73303E-7, true, true); //+/- 5sigma

dCutActionMap_BDTSignalCombo[locReactions[loc_i]]->Initialize(locEventLoop);

}

japp->RootWriteLock();

{

++Get_NumEventWriterThreads();

}

japp->RootUnLock();

}

DEventWriterROOT::~DEventWriterROOT(void)

{

japp->RootWriteLock();

{

--Get_NumEventWriterThreads();

if(Get_NumEventWriterThreads() != 0)

{

japp->RootUnLock();

return;

}

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

{

Get_OutputROOTFiles()[loc_i]->Write(0, TObject::kOverwrite);

Get_OutputROOTFiles()[loc_i]->Close();

delete Get_OutputROOTFiles()[loc_i];

}

Get_OutputROOTFiles().clear();

}

japp->RootUnLock();

}

int& DEventWriterROOT::Get_NumEventWriterThreads(void) const

{

static int locNumEventWriterThreads = 0;

return locNumEventWriterThreads;

}

string& DEventWriterROOT::Get_ThrownTreeFileName(void) const

{

//static so that it's not a member: can be changed in a call to a const function //object is const when the user gets it

static string locThrownTreeFileName = "";

return locThrownTreeFileName;

}

map<TTree*, map<string, TClonesArray*> >& DEventWriterROOT::Get_ClonesArrayMap(void) const

{

static map<TTree*, map<string, TClonesArray*> > locClonesArrayMap;

return locClonesArrayMap;

}

map<TTree*, map<string, TObject*> >& DEventWriterROOT::Get_TObjectMap(void) const

{

static map<TTree*, map<string, TObject*> > locTObjectMap;

return locTObjectMap;

}

map<TTree*, map<string, unsigned int> >& DEventWriterROOT::Get_FundamentalArraySizeMap(void) const

{

static map<TTree*, map<string, unsigned int> > locFundamentalArraySizeMap;

return locFundamentalArraySizeMap;

}

deque<TFile*>& DEventWriterROOT::Get_OutputROOTFiles(void) const

{

static deque<TFile*> locOutputROOTFiles;

100

return locOutputROOTFiles;

101

}

102

103

void DEventWriterROOT::Create_ThrownTree(string locOutputFileName) const

104

{

105

japp->RootWriteLock();

106

{

107

if(Get_ThrownTreeFileName() != "")

108

{

109

japp->RootUnLock();

110

return; //already created

111

}

112

Get_ThrownTreeFileName() = locOutputFileName;

113

114

//create ROOT file if it doesn't exist already

115

TFile* locFile = (TFile*)gROOT->FindObject(locOutputFileName.c_str());

116

if(locFile == NULL__null)

117

{

118

locFile = new TFile(locOutputFileName.c_str(), "RECREATE");

119

Get_OutputROOTFiles().push_back(locFile);

120

}

121

122

//create tree if it doesn't exist already

123

string locTreeName = "Thrown_Tree";

124

locFile->cd();

125

if(gDirectory(TDirectory::CurrentDirectory())->Get(locTreeName.c_str()) != NULL__null)

126

{

127

japp->RootUnLock();

128

return; //already created by another thread

129

}

130

TTree* locTree = new TTree(locTreeName.c_str(), locTreeName.c_str());

131

132

/******************************************************************** Create Branches ********************************************************************/

133

134

//create basic/misc. tree branches (run#, event#, etc.)

135

Create_Branch_Fundamental<UInt_t>(locTree, "RunNumber");

136

Create_Branch_Fundamental<ULong64_t>(locTree, "EventNumber");

137

138

//Thrown Data

139

Create_Branches_Thrown(locTree, true);

140

141

//CUSTOM

142

Create_CustomBranches_ThrownTree(locTree);

143

}

144

japp->RootUnLock();

145

}

146

147

void DEventWriterROOT::Create_DataTrees(JEventLoop* locEventLoop) const

148

{

149

vector<const DMCThrown*> locMCThrowns;

150

locEventLoop->Get(locMCThrowns);

151

152

vector<const DReaction*> locReactions;

153

Get_Reactions(locEventLoop, locReactions);

154

155

//CREATE TTREES

156

japp->RootWriteLock();

157

{

158

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

159

{

160

if(locReactions[loc_i]->Get_EnableTTreeOutputFlag())

161

Create_DataTree(locReactions[loc_i], !locMCThrowns.empty());

162

}

163

}

164

japp->RootUnLock();

165

}

166

167

void DEventWriterROOT::Create_DataTree(const DReaction* locReaction, bool locIsMCDataFlag) const

168

{

169

string locReactionName = locReaction->Get_ReactionName();

170

171

//create ROOT file if it doesn't exist already

172

string locOutputFileName = locReaction->Get_TTreeOutputFileName();

173

TFile* locFile = (TFile*)gROOT->FindObject(locOutputFileName.c_str());

174

if(locFile == NULL__null)

175

{

176

locFile = new TFile(locOutputFileName.c_str(), "RECREATE");

177

Get_OutputROOTFiles().push_back(locFile);

178

}

179

180

//create tree if it doesn't exist already

181

string locTreeName = locReactionName + string("_Tree");

182

locFile->cd();

183

if(gDirectory(TDirectory::CurrentDirectory())->Get(locTreeName.c_str()) != NULL__null)

184

return; //already created by another thread

185

TTree* locTree = new TTree(locTreeName.c_str(), locTreeName.c_str());

186

187

//find the # particles of each pid

188

map<Particle_t, unsigned int> locParticleNumberMap;

189

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

190

{

191

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

192

deque<Particle_t> locFinalParticleIDs;

193

locReactionStep->Get_FinalParticleIDs(locFinalParticleIDs);

194

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

195

{

196

if(locReactionStep->Get_MissingParticleIndex() == int(loc_j)) //missing particle

197

continue;

198

Particle_t locPID = locFinalParticleIDs[loc_j];

199

if(locParticleNumberMap.find(locPID) == locParticleNumberMap.end())

200

locParticleNumberMap[locPID] = 1;

201

else

202

++locParticleNumberMap[locPID];

203

}

204

}

205

206

//fill maps

207

Create_UserInfoMaps(locTree, locReaction, locParticleNumberMap);

208

209

/******************************************************************** Create Branches ********************************************************************/

210

211

//create basic/misc. tree branches (run#, event#, etc.)

212

Create_Branch_Fundamental<UInt_t>(locTree, "RunNumber");

213

Create_Branch_Fundamental<ULong64_t>(locTree, "EventNumber");

214

215

//create thrown branches

216

if(locIsMCDataFlag)

217

{

218

Create_Branches_Thrown(locTree, false);

219

Create_Branch_Fundamental<Float_t>(locTree, "IsThrownTopology");

220

}

221

222

Particle_t locInitialPID = locReaction->Get_ReactionStep(0)->Get_InitialParticleID();

223

bool locBeamUsedFlag = ((locInitialPID == Gamma) || (locInitialPID == Electron) || (locInitialPID == Positron));

224

225

//create branches for final-state particle hypotheses

226

if(locBeamUsedFlag)

227

Create_Branches_Beam(locTree, locIsMCDataFlag);

228

Create_Branches_NeutralShowers(locTree, locIsMCDataFlag);

229

Create_Branches_ChargedHypotheses(locTree, locIsMCDataFlag);

230

231

//create branches for combos

232

Create_Branches_Combo(locTree, locReaction, locIsMCDataFlag, locParticleNumberMap);

233

234

//Custom branches

235

Create_CustomBranches_DataTree(locTree, locReaction, locIsMCDataFlag);

236

}

237

238

void DEventWriterROOT::Create_UserInfoMaps(TTree* locTree, const DReaction* locReaction, map<Particle_t, unsigned int>& locParticleNumberMap) const

239

{

240

//kinfit type

241

DKinFitType locKinFitType = locReaction->Get_KinFitType();

242

243

//create & add reaction identification maps

244

TList* locUserInfo = locTree->GetUserInfo();

245

TMap* locNameToPIDMap = new TMap();

246

locNameToPIDMap->SetName("NameToPIDMap");

247

locUserInfo->Add(locNameToPIDMap);

248

249

TMap* locNameToPositionMap = new TMap(); //not filled for target or initial particles

250

locNameToPositionMap->SetName("NameToPositionMap");

251

locUserInfo->Add(locNameToPositionMap);

252

253

TMap* locPositionToNameMap = new TMap(); //not filled for target or initial particles

254

locPositionToNameMap->SetName("PositionToNameMap");

255

locUserInfo->Add(locPositionToNameMap);

256

257

TMap* locPositionToPIDMap = new TMap();

258

locPositionToPIDMap->SetName("PositionToPIDMap");

259

locUserInfo->Add(locPositionToPIDMap);

260

261

TMap* locDecayProductMap = new TMap(); //excludes resonances!!! //excludes intermediate decays (e.g. if xi- -> pi-, lambda -> pi-, pi-, p: will be xi- -> pi-, pi-, p and no lambda decay present)

262

locDecayProductMap->SetName("DecayProductMap"); //parent name string -> tlist of decay product name strings

263

locUserInfo->Add(locDecayProductMap);

264

265

TMap* locMiscInfoMap = new TMap();

266

locMiscInfoMap->SetName("MiscInfoMap");

267

locUserInfo->Add(locMiscInfoMap);

268

269

TList* locParticleNameList = new TList();

270

locParticleNameList->SetName("ParticleNameList");

271

locUserInfo->Add(locParticleNameList);

272

273

ostringstream locKinFitTypeStream;

274

locKinFitTypeStream << locKinFitType;

275

locMiscInfoMap->Add(new TObjString("KinFitType"), new TObjString(locKinFitTypeStream.str().c_str()));

276

277

map<Particle_t, unsigned int> locParticleNumberMap_Current;

278

Particle_t locPID;

279

TObjString *locObjString_PID, *locObjString_Position, *locObjString_ParticleName;

280

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

281

{

282

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

283

284

//initial particle

285

{

286

ostringstream locPIDStream, locPositionStream, locParticleNameStream;

287

locPID = locReactionStep->Get_InitialParticleID();

288

locPIDStream << PDGtype(locPID);

289

locObjString_PID = new TObjString(locPIDStream.str().c_str());

290

locPositionStream << loc_i << "_" << -1;

291

locObjString_Position = new TObjString(locPositionStream.str().c_str());

292

locPositionToPIDMap->Add(locObjString_Position, locObjString_PID);

293

if((loc_i == 0) && ((locPID == Gamma) || (locPID == Electron) || (locPID == Positron)))

294

{

295

locParticleNameStream << "ComboBeam";

296

locObjString_ParticleName = new TObjString(locParticleNameStream.str().c_str());

297

locNameToPIDMap->Add(locObjString_ParticleName, locObjString_PID);

298

locParticleNameList->AddLast(locObjString_ParticleName);

299

}

300

}

301

302

//target particle

303

locPID = locReactionStep->Get_TargetParticleID();

304

if((loc_i == 0) && (locPID != Unknown))

305

{

306

ostringstream locPIDStream, locPositionStream, locParticleNameStream;

307

locPIDStream << PDGtype(locPID);

308

locObjString_PID = new TObjString(locPIDStream.str().c_str());

309

locPositionStream << loc_i << "_" << -2;

310

locObjString_Position = new TObjString(locPositionStream.str().c_str());

311

locPositionToPIDMap->Add(locObjString_Position, locObjString_PID);

312

locParticleNameStream << "Target" << Convert_ToBranchName(ParticleType(locPID));

313

locObjString_ParticleName = new TObjString(locParticleNameStream.str().c_str());

314

locNameToPositionMap->Add(locObjString_ParticleName, locObjString_Position);

315

locNameToPIDMap->Add(locObjString_ParticleName, locObjString_PID);

316

string locPIDName = locParticleNameStream.str() + string("__PID");

317

locMiscInfoMap->Add(new TObjString(locPIDName.c_str()), locObjString_PID);

318

ostringstream locMassStream;

319

locMassStream << ParticleMass(locPID);

320

locMiscInfoMap->Add(new TObjString("Target__Mass"), new TObjString(locMassStream.str().c_str()));

321

locParticleNameList->AddLast(locObjString_ParticleName);

322

}

323

324

//final particles

325

deque<Particle_t> locFinalParticleIDs;

326

locReactionStep->Get_FinalParticleIDs(locFinalParticleIDs);

327

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

328

{

329

ostringstream locPIDStream, locPositionStream;

330

locPID = locFinalParticleIDs[loc_j];

331

locPIDStream << PDGtype(locPID);

332

locObjString_PID = new TObjString(locPIDStream.str().c_str());

333

334

locPositionStream << loc_i << "_" << loc_j;

335

locObjString_Position = new TObjString(locPositionStream.str().c_str());

336

337

if(locReactionStep->Get_MissingParticleIndex() == int(loc_j)) //missing particle

338

{

339

ostringstream locParticleNameStream;

340

locParticleNameStream << "Missing" << Convert_ToBranchName(ParticleType(locPID));

341

locObjString_ParticleName = new TObjString(locParticleNameStream.str().c_str());

342

locNameToPositionMap->Add(locObjString_ParticleName, locObjString_Position);

343

locPositionToNameMap->Add(locObjString_Position, locObjString_ParticleName);

344

locNameToPIDMap->Add(locObjString_ParticleName, locObjString_PID);

345

string locPIDName = locParticleNameStream.str() + string("__PID");

346

locMiscInfoMap->Add(new TObjString(locPIDName.c_str()), locObjString_PID);

347

ostringstream locMassStream;

348

locMassStream << ParticleMass(locPID);

349

string locMassName = locParticleNameStream.str() + string("__Mass");

350

locMiscInfoMap->Add(new TObjString(locMassName.c_str()), new TObjString(locMassStream.str().c_str()));

351

locParticleNameList->AddLast(locObjString_ParticleName);

352

continue;

353

}

354

355

if(locParticleNumberMap_Current.find(locPID) == locParticleNumberMap_Current.end())

356

locParticleNumberMap_Current[locPID] = 1;

357

else

358

++locParticleNumberMap_Current[locPID];

359

360

//see if decays in a future step

361

bool locDecaysFlag = false;

362

for(size_t loc_k = loc_i + 1; loc_k < locReaction->Get_NumReactionSteps(); ++loc_k)

363

{

364

if(locReaction->Get_ReactionStep(loc_k)->Get_InitialParticleID() != locPID)

365

continue;

366

locDecaysFlag = true;

367

break;

368

}

369

370

ostringstream locParticleNameStream;

371

if(locDecaysFlag)

372

locParticleNameStream << "Decaying";

373

locParticleNameStream << Convert_ToBranchName(ParticleType(locPID));

374

if(locParticleNumberMap[locPID] > 1)

375

locParticleNameStream << locParticleNumberMap_Current[locPID];

376

locObjString_ParticleName = new TObjString(locParticleNameStream.str().c_str());

377

locParticleNameList->AddLast(locObjString_ParticleName);

378

379

locPositionToPIDMap->Add(locObjString_Position, locObjString_PID);

380

locNameToPositionMap->Add(locObjString_ParticleName, locObjString_Position);

381

locPositionToNameMap->Add(locObjString_Position, locObjString_ParticleName);

382

locNameToPIDMap->Add(locObjString_ParticleName, locObjString_PID);

383

if(locDecaysFlag)

384

{

385

ostringstream locMassStream;

386

locMassStream << ParticleMass(locPID);

387

string locMassName = locParticleNameStream.str() + string("__Mass");

388

locMiscInfoMap->Add(new TObjString(locMassName.c_str()), new TObjString(locMassStream.str().c_str()));

389

}

390

}

391

}

392

393

//fill decay product map

394

deque<size_t> locSavedSteps;

395

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

396

{

397

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

398

399

//initial particle

400

locPID = locReactionStep->Get_InitialParticleID();

401

if(loc_i == 0)

402

continue;

403

if(!IsFixedMass(locPID))

404

continue; //don't save resonance decays to the map

405

406

//check to see if this decay has already been saved

407

bool locStepAlreadySavedFlag = false;

408

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

409

{

410

if(locSavedSteps[loc_j] != loc_i)

411

continue;

412

locStepAlreadySavedFlag = true;

413

break;

414

}

415

if(locStepAlreadySavedFlag)

416

continue;

417

418

//find the name of the decay parent

419

//first need to find what instance this pid is a decay parent

420

size_t locDecayParentInstance = 1;

421

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

422

{

423

if(locReaction->Get_ReactionStep(loc_j)->Get_InitialParticleID() == locPID)

424

++locDecayParentInstance;

425

}

426

//construct the name

427

ostringstream locParticleNameStream;

428

locParticleNameStream << "Decaying";

429

locParticleNameStream << Convert_ToBranchName(ParticleType(locPID));

430

if(locParticleNumberMap[locPID] > 1)

431

locParticleNameStream << locDecayParentInstance;

432

locObjString_ParticleName = new TObjString(locParticleNameStream.str().c_str());

433

434

TList* locDecayProductNames = NULL__null;

435

Get_DecayProductNames(locReaction, loc_i, locPositionToNameMap, locDecayProductNames, locSavedSteps);

436

locDecayProductMap->Add(locObjString_ParticleName, locDecayProductNames); //parent name string -> tobjarray of decay product name strings

437

}

438

}

439

440

void DEventWriterROOT::Get_DecayProductNames(const DReaction* locReaction, size_t locReactionStepIndex, TMap* locPositionToNameMap, TList*& locDecayProductNames, deque<size_t>& locSavedSteps) const

441

{

442

const DReactionStep* locReactionStep = locReaction->Get_ReactionStep(locReactionStepIndex);

443

444

if(locDecayProductNames == NULL__null)

445

locDecayProductNames = new TList();

446

447

deque<Particle_t> locFinalParticleIDs;

448

locReactionStep->Get_FinalParticleIDs(locFinalParticleIDs);

449

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

450

{

451

Particle_t locPID = locFinalParticleIDs[loc_j];

452

453

//see if decays in a future step

454

bool locDecaysFlag = false;

455

for(size_t loc_k = locReactionStepIndex + 1; loc_k < locReaction->Get_NumReactionSteps(); ++loc_k)

456

{

457

if(locReaction->Get_ReactionStep(loc_k)->Get_InitialParticleID() != locPID)

458

continue;

459

locDecaysFlag = true;

460

break;

461

}

462

463

//save and continue if doesn't decay

464

if(!locDecaysFlag)

465

{

466

ostringstream locPositionStream;

467

locPositionStream << locReactionStepIndex << "_" << loc_j;

468

locDecayProductNames->AddLast(locPositionToNameMap->GetValue(locPositionStream.str().c_str()));

469

continue;

470

}

471

472

//decays: find step where it decays at

473

//first find what instance in the final state this particle is

474

size_t locFinalStateInstance = 1;

475

for(size_t loc_k = 0; loc_k < locReactionStepIndex; ++loc_k)

476

{

477

const DReactionStep* locNewReactionStep = locReaction->Get_ReactionStep(loc_k);

478

deque<Particle_t> locNewFinalParticleIDs;

479

locNewReactionStep->Get_FinalParticleIDs(locNewFinalParticleIDs);

480

for(size_t loc_l = 0; loc_l < locNewFinalParticleIDs.size(); ++loc_l)

481

{

482

if((loc_k == locReactionStepIndex) && (loc_l == loc_j))

483

break;

484

if(locNewFinalParticleIDs[loc_l] == locPID)

485

++locFinalStateInstance;

486

}

487

}

488

489

//next find the step index that this particle decays at

490

size_t locNewReactionStepIndex = 0;

491

size_t locInitialStateInstance = 0;

492

for(size_t loc_k = 0; loc_k < locReaction->Get_NumReactionSteps(); ++loc_k)

493

{

494

if(locReaction->Get_ReactionStep(loc_k)->Get_InitialParticleID() != locPID)

495

continue;

496

++locInitialStateInstance;

497

if(locInitialStateInstance != locFinalStateInstance)

498

continue;

499

locNewReactionStepIndex = loc_k;

500

break;

501

}

502

503

//add decay products

504

Get_DecayProductNames(locReaction, locNewReactionStepIndex, locPositionToNameMap, locDecayProductNames, locSavedSteps);

505

}

506

507

locSavedSteps.push_back(locReactionStepIndex);

508

}

509

510

void DEventWriterROOT::Create_Branches_Thrown(TTree* locTree, bool locIsOnlyThrownFlag) const

511

{

512

//BEAM

513

Create_Branch_Fundamental<Int_t>(locTree, "ThrownBeam", "PID");

514

Create_Branch_NoSplitTObject<TLorentzVector>(locTree, "ThrownBeam", "X4"); //reported at target center

515

Create_Branch_NoSplitTObject<TLorentzVector>(locTree, "ThrownBeam", "P4");

516

517

//EVENT-WIDE INFO

518

Create_Branch_Fundamental<ULong64_t>(locTree, "NumPIDThrown_FinalState"); //19 digits

519

Create_Branch_Fundamental<ULong64_t>(locTree, "PIDThrown_Decaying");

520

Create_Branch_Fundamental<Float_t>(locTree, "MCWeight");

521

522

//PRODUCTS

523

Create_Branches_ThrownParticles(locTree, locIsOnlyThrownFlag);

524

}

525

526

void DEventWriterROOT::Create_Branches_ThrownParticles(TTree* locTree, bool locIsOnlyThrownFlag) const

527

{

528

string locParticleBranchName = "Thrown";

529

530

string locArraySizeString = "NumThrown";

531

Create_Branch_Fundamental<UInt_t>(locTree, locArraySizeString);

532

533

//IDENTIFIERS

534

Create_Branch_FundamentalArray<Int_t>(locTree, locParticleBranchName, "ParentIndex", locArraySizeString, dInitNumThrownArraySize);

535

Create_Branch_FundamentalArray<Int_t>(locTree, locParticleBranchName, "PID", locArraySizeString, dInitNumThrownArraySize);

536

if(!locIsOnlyThrownFlag)

537

{

538

Create_Branch_FundamentalArray<Int_t>(locTree, locParticleBranchName, "MatchID", locArraySizeString, dInitNumThrownArraySize);

539

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "MatchFOM", locArraySizeString, dInitNumThrownArraySize);

540

}

541

542

//KINEMATICS: THROWN //at the production vertex

543

Create_Branch_ClonesArray(locTree, locParticleBranchName, "X4", "TLorentzVector", dInitNumThrownArraySize);

544

Create_Branch_ClonesArray(locTree, locParticleBranchName, "P4", "TLorentzVector", dInitNumThrownArraySize);

545

}

546

547

void DEventWriterROOT::Create_Branches_Beam(TTree* locTree, bool locIsMCDataFlag) const

548

{

549

string locArraySizeString = "NumBeam";

550

Create_Branch_Fundamental<UInt_t>(locTree, locArraySizeString);

551

552

//IDENTIFIER

553

Create_Branch_FundamentalArray<Int_t>(locTree, "Beam", "PID", locArraySizeString, dInitNumBeamArraySize);

554

if(locIsMCDataFlag)

555

Create_Branch_FundamentalArray<Bool_t>(locTree, "Beam", "IsGenerator", locArraySizeString, dInitNumBeamArraySize);

556

557

//KINEMATICS: MEASURED //at the target center

558

Create_Branch_ClonesArray(locTree, "Beam", "X4_TargetCenter", "TLorentzVector", dInitNumBeamArraySize);

559

Create_Branch_ClonesArray(locTree, "Beam", "P4_Measured", "TLorentzVector", dInitNumBeamArraySize);

560

}

561

562

void DEventWriterROOT::Create_Branches_ChargedHypotheses(TTree* locTree, bool locIsMCDataFlag) const

563

{

564

string locArraySizeString = "NumChargedHypos";

565

Create_Branch_Fundamental<UInt_t>(locTree, locArraySizeString);

566

567

string locParticleBranchName = "ChargedHypo";

568

569

//IDENTIFIERS / MATCHING

570

Create_Branch_FundamentalArray<Int_t>(locTree, locParticleBranchName, "TrackID", locArraySizeString, dInitNumTrackArraySize);

571

Create_Branch_FundamentalArray<Int_t>(locTree, locParticleBranchName, "PID", locArraySizeString, dInitNumTrackArraySize);

572

if(locIsMCDataFlag)

573

Create_Branch_FundamentalArray<Int_t>(locTree, locParticleBranchName, "ThrownIndex", locArraySizeString, dInitNumTrackArraySize);

574

575

//KINEMATICS //at the production vertex

576

Create_Branch_ClonesArray(locTree, locParticleBranchName, "X4_Measured", "TLorentzVector", dInitNumTrackArraySize);

577

Create_Branch_ClonesArray(locTree, locParticleBranchName, "P4_Measured", "TLorentzVector", dInitNumTrackArraySize);

578

579

//TRACKING INFO

580

Create_Branch_FundamentalArray<UInt_t>(locTree, locParticleBranchName, "NDF_Tracking", locArraySizeString, dInitNumTrackArraySize);

581

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "ChiSq_Tracking", locArraySizeString, dInitNumTrackArraySize);

582

Create_Branch_FundamentalArray<UInt_t>(locTree, locParticleBranchName, "NDF_DCdEdx", locArraySizeString, dInitNumTrackArraySize);

583

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "ChiSq_DCdEdx", locArraySizeString, dInitNumTrackArraySize);

584

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "dEdx_CDC", locArraySizeString, dInitNumTrackArraySize);

585

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "dEdx_FDC", locArraySizeString, dInitNumTrackArraySize);

586

587

//TIMING INFO

588

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "HitTime", locArraySizeString, dInitNumTrackArraySize);

589

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "RFDeltaTVar", locArraySizeString, dInitNumTrackArraySize);

590

591

//HIT ENERGY

592

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "dEdx_TOF", locArraySizeString, dInitNumTrackArraySize);

593

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "dEdx_ST", locArraySizeString, dInitNumTrackArraySize);

594

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "Energy_BCAL", locArraySizeString, dInitNumTrackArraySize);

595

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "Energy_FCAL", locArraySizeString, dInitNumTrackArraySize);

596

597

//SHOWER MATCHING:

598

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "TrackBCAL_DeltaPhi", locArraySizeString, dInitNumTrackArraySize);

599

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "TrackBCAL_DeltaZ", locArraySizeString, dInitNumTrackArraySize);

600

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "TrackFCAL_DOCA", locArraySizeString, dInitNumTrackArraySize);

601

}

602

603

void DEventWriterROOT::Create_Branches_NeutralShowers(TTree* locTree, bool locIsMCDataFlag) const

604

{

605

string locArraySizeString = "NumNeutralShowers";

606

string locParticleBranchName = "NeutralShower";

607

Create_Branch_Fundamental<UInt_t>(locTree, locArraySizeString);

608

609

//MATCHING

610

if(locIsMCDataFlag)

611

Create_Branch_FundamentalArray<Int_t>(locTree, locParticleBranchName, "ThrownIndex", locArraySizeString, dInitNumShowerArraySize);

612

613

//SHOWER INFO

614

Create_Branch_ClonesArray(locTree, locParticleBranchName, "X4_Shower", "TLorentzVector", dInitNumShowerArraySize);

615

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "Energy_BCAL", locArraySizeString, dInitNumShowerArraySize);

616

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "Energy_FCAL", locArraySizeString, dInitNumShowerArraySize);

617

618

//NEARBY TRACKS

619

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "TrackBCAL_DeltaPhi", locArraySizeString, dInitNumShowerArraySize);

620

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "TrackBCAL_DeltaZ", locArraySizeString, dInitNumShowerArraySize);

621

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "TrackFCAL_DOCA", locArraySizeString, dInitNumShowerArraySize);

622

623

//PHOTON PID INFO

624

//Computed using DVertex (best estimate of reaction vertex using all "good" tracks)

625

//Can be used to compute timing chisq //is invalid for non-photons, so computed assuming photon PID

626

//Variance of X4_Measured.T() - RFTime, regardless of which RF bunch is chosen. //RF bunch is combo-depende

627

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "PhotonRFDeltaTVar", locArraySizeString, dInitNumShowerArraySize);

628

}

629

630

void DEventWriterROOT::Create_Branches_Combo(TTree* locTree, const DReaction* locReaction, bool locIsMCDataFlag, const map<Particle_t, unsigned int>& locParticleNumberMap) const

631

{

632

string locNumComboString = "NumCombos";

633

Create_Branch_Fundamental<UInt_t>(locTree, locNumComboString);

634

635

//kinfit type

636

DKinFitType locKinFitType = locReaction->Get_KinFitType();

637

bool locKinFitFlag = (locKinFitType != d_NoFit);

638

639

//Is-cut

640

Create_Branch_FundamentalArray<Bool_t>(locTree, "IsComboCut", locNumComboString, dInitNumComboArraySize);

641

642

//create combo-dependent, particle-independent branches

643

if(locIsMCDataFlag)

644

{

645

Create_Branch_FundamentalArray<Bool_t>(locTree, "IsTrueCombo", locNumComboString, dInitNumComboArraySize);

646

Create_Branch_FundamentalArray<Bool_t>(locTree, "IsBDTSignalCombo", locNumComboString, dInitNumComboArraySize);

647

}

648

649

Create_Branch_FundamentalArray<Float_t>(locTree, "RFTime_Measured", locNumComboString, dInitNumComboArraySize);

650

if(locKinFitFlag)

651

{

652

Create_Branch_FundamentalArray<Float_t>(locTree, "ChiSq_KinFit", locNumComboString, dInitNumComboArraySize);

653

Create_Branch_FundamentalArray<UInt_t>(locTree, "NDF_KinFit", locNumComboString, dInitNumComboArraySize);

654

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

655

Create_Branch_FundamentalArray<Float_t>(locTree, "RFTime_KinFit", locNumComboString, dInitNumComboArraySize);

656

}

657

658

map<Particle_t, unsigned int> locParticleNumberMap_Current;

659

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

660

{

661

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

662

663

//initial particle

664

Particle_t locInitialPID = locReactionStep->Get_InitialParticleID();

665

//should check to make sure the beam particle isn't missing...

666

if((locInitialPID == Gamma) || (locInitialPID == Electron) || (locInitialPID == Positron))

667

Create_Branches_BeamComboParticle(locTree, locInitialPID, locKinFitType);

668

else //decaying

669

{

670

ostringstream locParticleNameStream;

671

locParticleNameStream << "Decaying" << Convert_ToBranchName(ParticleType(locInitialPID));

672

if(IsFixedMass(locInitialPID) && ((locKinFitType == d_P4Fit) || (locKinFitType == d_P4AndVertexFit) || (locKinFitType == d_P4AndSpacetimeFit)))

673

Create_Branch_ClonesArray(locTree, locParticleNameStream.str(), "P4_KinFit", "TLorentzVector", dInitNumComboArraySize);

674

if(IsDetachedVertex(locInitialPID))

675

Create_Branch_ClonesArray(locTree, locParticleNameStream.str(), "X4", "TLorentzVector", dInitNumComboArraySize);

676

}

677

678

//final particles

679

deque<Particle_t> locFinalParticleIDs;

680

locReactionStep->Get_FinalParticleIDs(locFinalParticleIDs);

681

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

682

{

683

Particle_t locPID = locFinalParticleIDs[loc_j];

684

//decaying particle

685

if(locReaction->Check_IsDecayingParticle(locPID, loc_i + 1))

686

continue;

687

688

//missing particle

689

if(locReactionStep->Get_MissingParticleIndex() == int(loc_j))

690

{

691

// missing particle

692

ostringstream locParticleNameStream;

693

locParticleNameStream << "Missing" << Convert_ToBranchName(ParticleType(locPID));

694

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

695

Create_Branch_ClonesArray(locTree, locParticleNameStream.str(), "P4_KinFit", "TLorentzVector", dInitNumComboArraySize);

696

continue;

697

}

698

699

//detected

700

if(locParticleNumberMap_Current.find(locPID) == locParticleNumberMap_Current.end())

701

locParticleNumberMap_Current[locPID] = 1;

702

else

703

++locParticleNumberMap_Current[locPID];

704

705

ostringstream locParticleNameStream;

706

locParticleNameStream << Convert_ToBranchName(ParticleType(locPID));

707

if(locParticleNumberMap.find(locPID)->second > 1)

708

locParticleNameStream << locParticleNumberMap_Current[locPID];

709

710

if(ParticleCharge(locPID) == 0)

711

Create_Branches_ComboNeutral(locTree, locParticleNameStream.str(), locKinFitType);

712

else

713

Create_Branches_ComboTrack(locTree, locParticleNameStream.str(), locKinFitType);

714

}

715

}

716

}

717

718

void DEventWriterROOT::Create_Branches_BeamComboParticle(TTree* locTree, Particle_t locBeamPID, DKinFitType locKinFitType) const

719

{

720

string locParticleBranchName = "ComboBeam";

721

string locArraySizeString = "NumCombos";

722

723

//IDENTIFIER

724

Create_Branch_FundamentalArray<Int_t>(locTree, locParticleBranchName, "BeamIndex", locArraySizeString, dInitNumComboArraySize);

725

726

//SPACETIME AT INTERACTION VERTEX

727

Create_Branch_ClonesArray(locTree, locParticleBranchName, "X4_Measured", "TLorentzVector", dInitNumComboArraySize);

728

729

//KINEMATICS: KINFIT //at the interaction vertex

730

if(locKinFitType != d_NoFit)

731

{

732

if(((locKinFitType != d_VertexFit) && (locKinFitType != d_SpacetimeFit)) || (ParticleCharge(locBeamPID) != 0))

733

Create_Branch_ClonesArray(locTree, locParticleBranchName, "P4_KinFit", "TLorentzVector", dInitNumComboArraySize);

734

if(locKinFitType != d_P4Fit)

735

Create_Branch_ClonesArray(locTree, locParticleBranchName, "X4_KinFit", "TLorentzVector", dInitNumComboArraySize);

736

}

737

}

738

739

void DEventWriterROOT::Create_Branches_ComboTrack(TTree* locTree, string locParticleBranchName, DKinFitType locKinFitType) const

740

{

741

string locArraySizeString = "NumCombos";

742

743

//IDENTIFIER

744

Create_Branch_FundamentalArray<Int_t>(locTree, locParticleBranchName, "ChargedIndex", locArraySizeString, dInitNumComboArraySize); // is index to neutral/charged hypo array if PID is q0/q+-

745

746

//PID QUALITY

747

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "Beta_Timing_Measured", locArraySizeString, dInitNumComboArraySize);

748

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "ChiSq_Timing_Measured", locArraySizeString, dInitNumComboArraySize);

749

Create_Branch_FundamentalArray<UInt_t>(locTree, locParticleBranchName, "NDF_Timing", locArraySizeString, dInitNumComboArraySize);

750

751

//KINFIT INFO //at the production vertex

752

if(locKinFitType != d_NoFit)

753

{

754

Create_Branch_ClonesArray(locTree, locParticleBranchName, "P4_KinFit", "TLorentzVector", dInitNumComboArraySize);

755

if(locKinFitType != d_P4Fit)

756

{

757

Create_Branch_ClonesArray(locTree, locParticleBranchName, "X4_KinFit", "TLorentzVector", dInitNumComboArraySize);

758

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "Beta_Timing_KinFit", locArraySizeString, dInitNumComboArraySize);

759

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "ChiSq_Timing_KinFit", locArraySizeString, dInitNumComboArraySize);

760

}

761

}

762

}

763

764

void DEventWriterROOT::Create_Branches_ComboNeutral(TTree* locTree, string locParticleBranchName, DKinFitType locKinFitType) const

765

{

766

string locArraySizeString = "NumCombos";

767

768

//IDENTIFIER

769

Create_Branch_FundamentalArray<Int_t>(locTree, locParticleBranchName, "ShowerIndex", locArraySizeString, dInitNumComboArraySize); // is index to neutral/charged hypo array if PID is q0/q+-

770

771

//KINEMATICS //is combo-dependent: P4 defined by X4, X4 defined by other tracks

772

Create_Branch_ClonesArray(locTree, locParticleBranchName, "X4_Measured", "TLorentzVector", dInitNumComboArraySize);

773

Create_Branch_ClonesArray(locTree, locParticleBranchName, "P4_Measured", "TLorentzVector", dInitNumComboArraySize);

774

775

//PID QUALITY

776

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "Beta_Timing_Measured", locArraySizeString, dInitNumComboArraySize);

777

if(locParticleBranchName.substr(0, 6) == "Photon")

778

{

779

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "ChiSq_Timing_Measured", locArraySizeString, dInitNumComboArraySize);

780

Create_Branch_FundamentalArray<UInt_t>(locTree, locParticleBranchName, "NDF_Timing", locArraySizeString, dInitNumComboArraySize);

781

}

782

783

//KINFIT INFO //at the production vertex

784

if(locKinFitType != d_NoFit)

785

{

786

Create_Branch_ClonesArray(locTree, locParticleBranchName, "P4_KinFit", "TLorentzVector", dInitNumComboArraySize);

787

if(locKinFitType != d_P4Fit)

788

{

789

Create_Branch_ClonesArray(locTree, locParticleBranchName, "X4_KinFit", "TLorentzVector", dInitNumComboArraySize);

790

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "Beta_Timing_KinFit", locArraySizeString, dInitNumComboArraySize);

791

if(locParticleBranchName.substr(0, 6) == "Photon")

792

Create_Branch_FundamentalArray<Float_t>(locTree, locParticleBranchName, "ChiSq_Timing_KinFit", locArraySizeString, dInitNumComboArraySize);

793

}

794

}

795

}

796

797

void DEventWriterROOT::Get_Reactions(jana::JEventLoop* locEventLoop, vector<const DReaction*>& locReactions) const

798

{

799

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

800

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

801

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

802

// be used at a time.)

803

locReactions.clear();

804

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

805

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

806

{

807

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

808

if(locFactory == NULL__null)

809

continue;

810

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

811

continue;

812

813

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

814

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

815

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

816

// overall list.

817

vector<const DReaction*> locReactionsSubset;

818

locFactory->Get(locReactionsSubset);

819

locReactions.insert(locReactions.end(), locReactionsSubset.begin(), locReactionsSubset.end());

820

}

821

}

822

823

void DEventWriterROOT::Fill_ThrownTree(JEventLoop* locEventLoop) const

824

{

825

vector<const DMCThrown*> locMCThrowns_FinalState;

826

locEventLoop->Get(locMCThrowns_FinalState, "FinalState");

827

828

vector<const DMCThrown*> locMCThrowns_Decaying;

829

locEventLoop->Get(locMCThrowns_Decaying, "Decaying");

830

831

const DMCReaction* locMCReaction = NULL__null;

832

locEventLoop->GetSingle(locMCReaction);

833

834

//Pre-compute before entering lock

835

ULong64_t locNumPIDThrown_FinalState = 0, locPIDThrown_Decaying = 0;

836

Compute_ThrownPIDInfo(locMCThrowns_FinalState, locMCThrowns_Decaying, locNumPIDThrown_FinalState, locPIDThrown_Decaying);

837

838

//Pre-compute before entering lock

839

vector<const DMCThrown*> locMCThrownsToSave;

840

map<const DMCThrown*, unsigned int> locThrownIndexMap;

841

Group_ThrownParticles(locMCThrowns_FinalState, locMCThrowns_Decaying, locMCThrownsToSave, locThrownIndexMap);

842

843

string locTreeName = "Thrown_Tree";

844

japp->RootWriteLock();

845

{

846

if(Get_ThrownTreeFileName() == "")

847

{

848

japp->RootUnLock();

849

return;

850

}

851

852

//get the tree

853

TFile* locFile = (TFile*)gROOT->FindObject(Get_ThrownTreeFileName().c_str());

854

locFile->cd("/");

855

TTree* locTree = (TTree*)gDirectory(TDirectory::CurrentDirectory())->Get(locTreeName.c_str());

856

if(locTree == NULL__null)

857

{

858

cout << "ERROR: OUTPUT ROOT TREE NOT CREATED (in DEventWriterROOT::Fill_ThrownTree()). SKIP FILLING. " << endl;

859

japp->RootUnLock();

860

return;

861

}

862

863

//clear the tclonesarry's

864

map<string, TClonesArray*>& locTreeMap_ClonesArray = Get_ClonesArrayMap()[locTree];

865

map<string, TClonesArray*>::iterator locClonesArrayMapIterator;

866

for(locClonesArrayMapIterator = locTreeMap_ClonesArray.begin(); locClonesArrayMapIterator != locTreeMap_ClonesArray.end(); ++locClonesArrayMapIterator)

867

locClonesArrayMapIterator->second->Clear();

868

869

//primary event info

870

Fill_FundamentalData<UInt_t>(locTree, "RunNumber", locEventLoop->GetJEvent().GetRunNumber());

871

Fill_FundamentalData<ULong64_t>(locTree, "EventNumber", locEventLoop->GetJEvent().GetEventNumber());

872

873

//throwns

874

Fill_ThrownInfo(locTree, locMCReaction, locMCThrownsToSave, locThrownIndexMap, locNumPIDThrown_FinalState, locPIDThrown_Decaying);

875

876

//Custom Branches

877

Fill_CustomBranches_ThrownTree(locTree, locMCReaction, locMCThrownsToSave);

878

879

locTree->Fill();

880

}

881

japp->RootUnLock();

882

}

883

884

void DEventWriterROOT::Fill_DataTrees(JEventLoop* locEventLoop, string locDReactionTag) const

885

{

886

if(locDReactionTag == "Thrown")

887

{

888

cout << "WARNING: CANNOT FILL THROWN TREE WITH THIS FUNCTION." << endl;

889

return;

890

}

891

892

vector<const DAnalysisResults*> locAnalysisResultsVector;

893

locEventLoop->Get(locAnalysisResultsVector);

894

895

vector<const DReaction*> locReactionsWithTag;

896

locEventLoop->Get(locReactionsWithTag, locDReactionTag.c_str());

897

898

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

899

{

900

deque<const DParticleCombo*> locPassedParticleCombos;

901

locAnalysisResultsVector[loc_i]->Get_PassedParticleCombos(locPassedParticleCombos);

902

if(locPassedParticleCombos.empty())

903

continue;

904

const DReaction* locReaction = locAnalysisResultsVector[loc_i]->Get_Reaction();

905

if(!locReaction->Get_EnableTTreeOutputFlag())

906

continue;

907

bool locReactionFoundFlag = false;

908

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

909

{

910

if(locReactionsWithTag[loc_j] != locReaction)

911

continue;

912

locReactionFoundFlag = true;

913

break;

914

}

915

if(!locReactionFoundFlag)

916

continue; //reaction not from this factory, continue

917

Fill_DataTree(locEventLoop, locReaction, locPassedParticleCombos);

918

}

919

}

920

921

void DEventWriterROOT::Fill_DataTree(JEventLoop* locEventLoop, const DReaction* locReaction, deque<const DParticleCombo*>& locParticleCombos) const

922

{

923

if(locReaction->Get_ReactionName() == "Thrown")

924

{

925

cout << "WARNING: CANNOT FILL THROWN TREE WITH THIS FUNCTION." << endl;

926

return;

927

}

928

929

if(!locReaction->Get_EnableTTreeOutputFlag())

930

{

931

cout << "WARNING: ROOT TTREE OUTPUT NOT ENABLED FOR THIS DREACTION (" << locReaction->Get_ReactionName() << ")" << endl;

932

return;

933

}

934

935

//GET THROWN INFO

936

vector<const DMCThrown*> locMCThrowns_FinalState;

937

locEventLoop->Get(locMCThrowns_FinalState, "FinalState");

938

939

vector<const DMCThrown*> locMCThrowns_Decaying;

940

locEventLoop->Get(locMCThrowns_Decaying, "Decaying");

941

942

vector<const DMCThrownMatching*> locMCThrownMatchingVector;

943

locEventLoop->Get(locMCThrownMatchingVector);

944

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

945

946

vector<const DMCReaction*> locMCReactions;

947

locEventLoop->Get(locMCReactions);

948

const DMCReaction* locMCReaction = locMCReactions.empty() ? NULL__null : locMCReactions[0];

949

950

//Pre-compute before entering lock

951

ULong64_t locNumPIDThrown_FinalState = 0, locPIDThrown_Decaying = 0;

952

Compute_ThrownPIDInfo(locMCThrowns_FinalState, locMCThrowns_Decaying, locNumPIDThrown_FinalState, locPIDThrown_Decaying);

953

954

//Pre-compute before entering lock

955

vector<const DMCThrown*> locMCThrownsToSave;

956

map<const DMCThrown*, unsigned int> locThrownIndexMap;

957

Group_ThrownParticles(locMCThrowns_FinalState, locMCThrowns_Decaying, locMCThrownsToSave, locThrownIndexMap);

958

959

//GET DETECTOR MATCHES

960

const DDetectorMatches* locDetectorMatches = NULL__null;

961

locEventLoop->GetSingle(locDetectorMatches);

962

963

//GET NEUTRAL INFO

964

vector<const DNeutralParticle*> locNeutralParticles;

965

locEventLoop->Get(locNeutralParticles, dShowerSelectionTag.c_str());

966

967

//GET CHARGED TRACK INFO

968

vector<const DChargedTrack*> locChargedTracks;

969

locEventLoop->Get(locChargedTracks, dTrackSelectionTag.c_str());

970

971

vector<const DChargedTrackHypothesis*> locComboChargedTrackHypotheses;

972

locEventLoop->Get(locComboChargedTrackHypotheses, "Combo");

973

974

//Search for hypotheses not in the PreSelect factory (i.e. PID not in REST file)

975

map<const DChargedTrack*, vector<const DChargedTrackHypothesis*> > locNewHypothesesMap;

976

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

977

{

978

const DChargedTrackHypothesis* locOrigChargedTrackHypothesis = NULL__null;

979

locComboChargedTrackHypotheses[loc_i]->GetSingle(locOrigChargedTrackHypothesis);

980

if(locOrigChargedTrackHypothesis != NULL__null)

981

continue; //PID in REST file

982

983

const DChargedTrack* locOrigChargedTrack = NULL__null;

984

locComboChargedTrackHypotheses[loc_i]->GetSingle(locOrigChargedTrack);

985

986

locNewHypothesesMap[locOrigChargedTrack].push_back(locComboChargedTrackHypotheses[loc_i]);

987

}

988

989

//Build vector of combo-independent hypotheses to save

990

vector<const DChargedTrackHypothesis*> locIndependentChargedTrackHypotheses;

991

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

992

{

993

const vector<const DChargedTrackHypothesis*>& locTrackHypoVector = locChargedTracks[loc_i]->dChargedTrackHypotheses;

994

vector<const DChargedTrackHypothesis*>& locNewTrackHypoVector = locNewHypothesesMap[locChargedTracks[loc_i]];

995

996

locIndependentChargedTrackHypotheses.insert(locIndependentChargedTrackHypotheses.end(), locTrackHypoVector.begin(), locTrackHypoVector.end());

997

locIndependentChargedTrackHypotheses.insert(locIndependentChargedTrackHypotheses.end(), locNewTrackHypoVector.begin(), locNewTrackHypoVector.end());

998

}

999

1000

//Check whether beam is used in the combo

1001

Particle_t locInitialPID = locReaction->Get_ReactionStep(0)->Get_InitialParticleID();

1002

bool locBeamUsedFlag = ((locInitialPID == Gamma) || (locInitialPID == Electron) || (locInitialPID == Positron));

1003

1004

//GET BEAM PHOTONS

1005

//however, only fill with beam particles that are in the combos

1006

set<const DBeamPhoton*> locBeamPhotonSet;

1007

vector<const DBeamPhoton*> locBeamPhotons;

1008

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

1009

{

1010

const DParticleComboStep* locParticleComboStep = locParticleCombos[loc_j]->Get_ParticleComboStep(0);

1011

const DKinematicData* locKinematicData = locParticleComboStep->Get_InitialParticle_Measured();

1012

if(locKinematicData == NULL__null)

1013

continue;

1014

const DBeamPhoton* locBeamPhoton = dynamic_cast<const DBeamPhoton*>(locKinematicData);

1015

if(locBeamPhoton == NULL__null)

1016

continue;

1017

if(locBeamPhotonSet.find(locBeamPhoton) != locBeamPhotonSet.end())

1018

continue; //already included

1019

locBeamPhotonSet.insert(locBeamPhoton);

1020

locBeamPhotons.push_back(locBeamPhoton);

1021

}

1022

1023

//create map of particles to array index

1024

//used for pointing combo particles to the appropriate array index

1025

map<string, map<oid_t, int> > locObjectToArrayIndexMap;

1026

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

1027

locObjectToArrayIndexMap["DBeamPhoton"][locBeamPhotons[loc_i]->id] = loc_i;

1028

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

1029

locObjectToArrayIndexMap["DChargedTrackHypothesis"][locIndependentChargedTrackHypotheses[loc_i]->id] = loc_i;

1030

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

1031

locObjectToArrayIndexMap["DNeutralShower"][locNeutralParticles[loc_i]->dNeutralShower->dShowerID] = loc_i;

1032

1033

//EXECUTE ANALYSIS ACTIONS (Outside of ROOT Lock!)

1034

Bool_t locIsThrownTopologyFlag = kFALSE;

1035

vector<Bool_t> locIsTrueComboFlags;

1036

vector<Bool_t> locIsBDTSignalComboFlags;

1037

if(locMCReaction != NULL__null)

1038

{

1039

DCutAction_ThrownTopology* locThrownTopologyAction = dCutActionMap_ThrownTopology.find(locReaction)->second;

1040

locIsThrownTopologyFlag = (*locThrownTopologyAction)(locEventLoop, NULL__null); //combo not used/needed

1041

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

1042

{

1043

DCutAction_TrueCombo* locTrueComboAction = dCutActionMap_TrueCombo.find(locReaction)->second;

1044

locIsTrueComboFlags.push_back((*locTrueComboAction)(locEventLoop, locParticleCombos[loc_i]));

1045

1046

DCutAction_BDTSignalCombo* locBDTSignalComboAction = dCutActionMap_BDTSignalCombo.find(locReaction)->second;

1047

locIsBDTSignalComboFlags.push_back((*locBDTSignalComboAction)(locEventLoop, locParticleCombos[loc_i]));

1048

}

1049

}

1050

1051

string locOutputFileName = locReaction->Get_TTreeOutputFileName();

1052

string locTreeName = locReaction->Get_ReactionName() + string("_Tree");

1053

1054

japp->RootWriteLock();

1055

{

1056

TFile* locFile = (TFile*)gROOT->FindObject(locOutputFileName.c_str());

1057

if(locFile == NULL__null)

1058

{

1059

cout << "ERROR: OUTPUT ROOT TREE FILE NOT CREATED (in DEventWriterROOT::Fill_DataTree()). SKIP FILLING. " << endl;

1060

japp->RootUnLock();

1061

return;

1062

}

1063

1064

//get the tree info

1065

locFile->cd("/");

1066

TTree* locTree = (TTree*)gDirectory(TDirectory::CurrentDirectory())->Get(locTreeName.c_str());

1067

if(locTree == NULL__null)

1068

{

1069

cout << "ERROR: OUTPUT ROOT TREE NOT CREATED (in DEventWriterROOT::Fill_DataTree()). SKIP FILLING. " << endl;

1070

japp->RootUnLock();

1071

return;

1072

}

1073

1074

//clear the tclonesarry's

1075

map<string, TClonesArray*>& locTreeMap_ClonesArray = Get_ClonesArrayMap()[locTree];

1076

map<string, TClonesArray*>::iterator locClonesArrayMapIterator;

1077

for(locClonesArrayMapIterator = locTreeMap_ClonesArray.begin(); locClonesArrayMapIterator != locTreeMap_ClonesArray.end(); ++locClonesArrayMapIterator)

1078

locClonesArrayMapIterator->second->Clear();

1079

1080

//PRIMARY EVENT INFO

1081

Fill_FundamentalData<UInt_t>(locTree, "RunNumber", locEventLoop->GetJEvent().GetRunNumber());

1082

Fill_FundamentalData<ULong64_t>(locTree, "EventNumber", locEventLoop->GetJEvent().GetEventNumber());

1083

1084

//THROWN INFORMATION

1085

if(locMCReaction != NULL__null)

1086

{

1087

Fill_ThrownInfo(locTree, locMCReaction, locMCThrownsToSave, locThrownIndexMap, locNumPIDThrown_FinalState, locPIDThrown_Decaying, locMCThrownMatching, locObjectToArrayIndexMap);

1088

Fill_FundamentalData<Bool_t>(locTree, "IsThrownTopology", locIsThrownTopologyFlag);

1089

}

1090

1091

//INDEPENDENT BEAM PARTICLES

1092

if(locBeamUsedFlag)

1093

{

1094

//however, only fill with beam particles that are in the combos

1095

Fill_FundamentalData<UInt_t>(locTree, "NumBeam", locBeamPhotons.size());

1096

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

1097

Fill_BeamData(locTree, loc_i, locBeamPhotons[loc_i], locMCThrownMatching);

1098

}

1099

1100

//INDEPENDENT CHARGED TRACKS

1101

Fill_FundamentalData<UInt_t>(locTree, "NumChargedHypos", locIndependentChargedTrackHypotheses.size());

1102

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

1103

Fill_ChargedHypo(locTree, loc_i, locIndependentChargedTrackHypotheses[loc_i], locMCThrownMatching, locThrownIndexMap, locDetectorMatches);

1104

1105

//INDEPENDENT NEUTRAL SHOWERS

1106

Fill_FundamentalData<UInt_t>(locTree, "NumNeutralShowers", locNeutralParticles.size());

1107

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

1108

{

1109

const DNeutralParticleHypothesis* locPhotonHypothesis = locNeutralParticles[loc_i]->Get_Hypothesis(Gamma);

1110

Fill_NeutralShower(locTree, loc_i, locPhotonHypothesis, locMCThrownMatching, locThrownIndexMap, locDetectorMatches);

1111

}

1112

1113

//COMBOS

1114

Fill_FundamentalData<UInt_t>(locTree, "NumCombos", locParticleCombos.size());

1115

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

1116

{

1117

Fill_ComboData(locTree, locParticleCombos[loc_i], loc_i, locObjectToArrayIndexMap);

1118

if(locMCReaction != NULL__null)

1119

{

1120

Fill_FundamentalData<Bool_t>(locTree, "IsTrueCombo", locIsTrueComboFlags[loc_i], loc_i);

1121

Fill_FundamentalData<Bool_t>(locTree, "IsBDTSignalCombo", locIsBDTSignalComboFlags[loc_i], loc_i);

1122

}

1123

}

1124

1125

//CUSTOM

1126

Fill_CustomBranches_DataTree(locTree, locMCReaction, locMCThrownsToSave, locMCThrownMatching, locDetectorMatches, locBeamPhotons, locIndependentChargedTrackHypotheses, locNeutralParticles, locParticleCombos);

1127

1128

//FILL

1129

locTree->Fill();

1130

}

1131

japp->RootUnLock();

1132

}

1133

1134

ULong64_t DEventWriterROOT::Calc_ParticleMultiplexID(Particle_t locPID) const

1135

{

1136

int locPower = ParticleMultiplexPower(locPID);

1137

if(locPower == -1)

←

Taking true branch

→

1138

return 0;

←

Returning zero

→

1139

1140

int locIsFinalStateInt = Is_FinalStateParticle(locPID);

1141

if(locPID == Pi0)

1142

locIsFinalStateInt = 1;

1143

1144

if(locIsFinalStateInt == 1) //decimal

1145

{

1146

ULong64_t locParticleMultiplexID = 1;

1147

for(int loc_i = 0; loc_i < locPower; ++loc_i)

1148

locParticleMultiplexID *= ULong64_t(10);

1149

return locParticleMultiplexID;

1150

}

1151

//decaying: binary

1152

return (ULong64_t(1) << ULong64_t(locPower)); //bit-shift

1153

}

1154

1155

void DEventWriterROOT::Compute_ThrownPIDInfo(const vector<const DMCThrown*>& locMCThrowns_FinalState, const vector<const DMCThrown*>& locMCThrowns_Decaying, ULong64_t& locNumPIDThrown_FinalState, ULong64_t& locPIDThrown_Decaying) const

1156

{

1157

//THROWN PARTICLES BY PID

1158

locNumPIDThrown_FinalState = 0;

1159

for(size_t loc_i = 0; loc_i < locMCThrowns_FinalState.size(); ++loc_i) //final state

Loop condition is false. Execution continues on line 1168

→

1160

{

1161

Particle_t locPID = locMCThrowns_FinalState[loc_i]->PID();

1162

ULong64_t locPIDMultiplexID = Calc_ParticleMultiplexID(locPID);

1163

unsigned int locCurrentNumParticles = (locNumPIDThrown_FinalState / locPIDMultiplexID) % 10;

1164

if(locCurrentNumParticles != 9)

1165

locNumPIDThrown_FinalState += locPIDMultiplexID;

1166

}

1167

1168

locPIDThrown_Decaying = 0;

1169

for(size_t loc_i = 0; loc_i < locMCThrowns_Decaying.size(); ++loc_i) //decaying

←

Loop condition is true. Entering loop body

→

1170

{

1171

Particle_t locPID = locMCThrowns_Decaying[loc_i]->PID();

1172

ULong64_t locPIDMultiplexID = Calc_ParticleMultiplexID(locPID);

←

Calling 'DEventWriterROOT::Calc_ParticleMultiplexID'

→

←

Returning from 'DEventWriterROOT::Calc_ParticleMultiplexID'

→

←

'locPIDMultiplexID' initialized to 0

→

1173

if(locPID != Pi0)

←

Assuming 'locPID' is equal to Pi0

→

←

Taking false branch

→

1174

locPIDThrown_Decaying |= locPIDMultiplexID; //bit-wise or

1175

else //save pi0's as final state instead of decaying

1176

{

1177

unsigned int locCurrentNumParticles = (locNumPIDThrown_FinalState / locPIDMultiplexID) % 10;

←

Division by zero

1178

if(locCurrentNumParticles != 9)

1179

locNumPIDThrown_FinalState += locPIDMultiplexID;

1180

}

1181

}

1182

}

1183

1184

void DEventWriterROOT::Group_ThrownParticles(const vector<const DMCThrown*>& locMCThrowns_FinalState, const vector<const DMCThrown*>& locMCThrowns_Decaying, vector<const DMCThrown*>& locMCThrownsToSave, map<const DMCThrown*, unsigned int>& locThrownIndexMap) const

1185

{

1186

locMCThrownsToSave.clear();

1187

locMCThrownsToSave.insert(locMCThrownsToSave.end(), locMCThrowns_FinalState.begin(), locMCThrowns_FinalState.end());

1188

locMCThrownsToSave.insert(locMCThrownsToSave.end(), locMCThrowns_Decaying.begin(), locMCThrowns_Decaying.end());

1189

1190

//create map of mcthrown to array index

1191

locThrownIndexMap.clear();

1192

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

1193

locThrownIndexMap[locMCThrownsToSave[loc_i]] = loc_i;

1194

}

1195

1196

void DEventWriterROOT::Fill_ThrownInfo(TTree* locTree, const DMCReaction* locMCReaction, const vector<const DMCThrown*>& locMCThrowns, const map<const DMCThrown*, unsigned int>& locThrownIndexMap, ULong64_t locNumPIDThrown_FinalState, ULong64_t locPIDThrown_Decaying) const

1197

{

1198

map<string, map<oid_t, int> > locObjectToArrayIndexMap;

1199

Fill_ThrownInfo(locTree, locMCReaction, locMCThrowns, locThrownIndexMap, locNumPIDThrown_FinalState, locPIDThrown_Decaying, NULL__null, locObjectToArrayIndexMap);

1200

}

1201

1202

1203

{

1204

//THIS MUST BE CALLED FROM WITHIN A LOCK, SO DO NOT PASS IN JEVENTLOOP! //TOO TEMPTING TO DO SOMETHING BAD

1205

1206

//WEIGHT

1207

Fill_FundamentalData<Float_t>(locTree, "MCWeight", locMCReaction->weight);

1208

1209

//THROWN BEAM

1210

Fill_FundamentalData<Int_t>(locTree, "ThrownBeam", "PID", PDGtype(locMCReaction->beam.PID()));

1211

1212

DVector3 locThrownBeamX3 = locMCReaction->beam.position();

1213

TLorentzVector locThrownBeamTX4(locThrownBeamX3.X(), locThrownBeamX3.Y(), locThrownBeamX3.Z(), locMCReaction->beam.time());

1214

Fill_TObjectData<TLorentzVector>(locTree, "ThrownBeam", "X4", locThrownBeamTX4);

1215

1216

DLorentzVector locThrownBeamP4 = locMCReaction->beam.lorentzMomentum();

1217

TLorentzVector locThrownBeamTP4(locThrownBeamP4.Px(), locThrownBeamP4.Py(), locThrownBeamP4.Pz(), locThrownBeamP4.E());

1218

Fill_TObjectData<TLorentzVector>(locTree, "ThrownBeam", "P4", locThrownBeamTP4);

1219

1220

//THROWN PRODUCTS

1221

Fill_FundamentalData<UInt_t>(locTree, "NumThrown", locMCThrowns.size());

1222

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

1223

Fill_ThrownParticleData(locTree, loc_i, locMCThrowns[loc_i], locThrownIndexMap, locMCThrownMatching, locObjectToArrayIndexMap);

1224

1225

//PID INFO

1226

Fill_FundamentalData<ULong64_t>(locTree, "NumPIDThrown_FinalState", locNumPIDThrown_FinalState); //19 digits

1227

Fill_FundamentalData<ULong64_t>(locTree, "PIDThrown_Decaying", locPIDThrown_Decaying);

1228

}

1229

1230

void DEventWriterROOT::Fill_ThrownParticleData(TTree* locTree, unsigned int locArrayIndex, const DMCThrown* locMCThrown, const map<const DMCThrown*, unsigned int>& locThrownIndexMap, const DMCThrownMatching* locMCThrownMatching, const map<string, map<oid_t, int> >& locObjectToArrayIndexMap) const

1231

{

1232

string locParticleBranchName = "Thrown";

1233

1234

//IDENTIFIERS

1235

int locParentIndex = -1; //e.g. photoproduced

1236

map<const DMCThrown*, unsigned int>::const_iterator locIterator;

1237

for(locIterator = locThrownIndexMap.begin(); locIterator != locThrownIndexMap.end(); ++locIterator)

1238

{

1239

if(locIterator->first->myid != locMCThrown->parentid)

1240

continue;

1241

locParentIndex = locIterator->second;

1242

break;

1243

}

1244

Fill_FundamentalData<Int_t>(locTree, locParticleBranchName, "ParentIndex", locParentIndex, locArrayIndex);

1245

Fill_FundamentalData<Int_t>(locTree, locParticleBranchName, "PID", locMCThrown->pdgtype, locArrayIndex);

1246

1247

//MATCHING

1248

if(locMCThrownMatching != NULL__null)

1249

{

1250

Int_t locMatchID = -1;

1251

double locMatchFOM = -1.0;

1252

if(ParticleCharge(locMCThrown->PID()) != 0)

1253

{

1254

const DChargedTrack* locChargedTrack = locMCThrownMatching->Get_MatchingChargedTrack(locMCThrown, locMatchFOM);

1255

if(locChargedTrack != NULL__null)

1256

locMatchID = locChargedTrack->Get_BestFOM()->candidateid;

1257

}

1258

else

1259

{

1260

//Can't use DNeutralShower JObject::id:

1261

//Matching done with default-tag showers, but pre-select showers are saved to tree: JObject::id's don't match

1262

const DNeutralShower* locNeutralShower = locMCThrownMatching->Get_MatchingNeutralShower(locMCThrown, locMatchFOM);

1263

if(locNeutralShower != NULL__null)

1264

locMatchID = locObjectToArrayIndexMap.find("DNeutralShower")->second.find(locNeutralShower->dShowerID)->second;

1265

}

1266

Fill_FundamentalData<Int_t>(locTree, locParticleBranchName, "MatchID", locMatchID, locArrayIndex);

1267

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "MatchFOM", locMatchFOM, locArrayIndex);

1268

}

1269

1270

//KINEMATICS: THROWN //at the production vertex

1271

TLorentzVector locX4_Thrown(locMCThrown->position().X(), locMCThrown->position().Y(), locMCThrown->position().Z(), locMCThrown->time());

1272

Fill_ClonesData<TLorentzVector>(locTree, locParticleBranchName, "X4", locX4_Thrown, locArrayIndex);

1273

TLorentzVector locP4_Thrown(locMCThrown->momentum().X(), locMCThrown->momentum().Y(), locMCThrown->momentum().Z(), locMCThrown->energy());

1274

Fill_ClonesData<TLorentzVector>(locTree, locParticleBranchName, "P4", locP4_Thrown, locArrayIndex);

1275

}

1276

1277

void DEventWriterROOT::Fill_BeamData(TTree* locTree, unsigned int locArrayIndex, const DBeamPhoton* locBeamPhoton, const DMCThrownMatching* locMCThrownMatching) const

1278

{

1279

string locParticleBranchName = "Beam";

1280

1281

//IDENTIFIER

1282

Fill_FundamentalData<Int_t>(locTree, locParticleBranchName, "PID", PDGtype(locBeamPhoton->PID()), locArrayIndex);

1283

1284

//MATCHING

1285

if(locMCThrownMatching != NULL__null)

1286

{

1287

Bool_t locIsGeneratorFlag = (locMCThrownMatching->Get_ReconMCGENBeamPhoton() == locBeamPhoton) ? kTRUE : kFALSE;

1288

Fill_FundamentalData<Bool_t>(locTree, locParticleBranchName, "IsGenerator", locIsGeneratorFlag, locArrayIndex);

1289

}

1290

1291

//KINEMATICS: MEASURED

1292

DVector3 locPosition = locBeamPhoton->position();

1293

TLorentzVector locX4_Measured(locPosition.X(), locPosition.Y(), locPosition.Z(), locBeamPhoton->time());

1294

Fill_ClonesData<TLorentzVector>(locTree, locParticleBranchName, "X4_TargetCenter", locX4_Measured, locArrayIndex);

1295

1296

DLorentzVector locDP4 = locBeamPhoton->lorentzMomentum();

1297

TLorentzVector locP4_Measured(locDP4.Px(), locDP4.Py(), locDP4.Pz(), locDP4.E());

1298

Fill_ClonesData<TLorentzVector>(locTree, locParticleBranchName, "P4_Measured", locP4_Measured, locArrayIndex);

1299

}

1300

1301

void DEventWriterROOT::Fill_ChargedHypo(TTree* locTree, unsigned int locArrayIndex, const DChargedTrackHypothesis* locChargedTrackHypothesis, const DMCThrownMatching* locMCThrownMatching, const map<const DMCThrown*, unsigned int>& locThrownIndexMap, const DDetectorMatches* locDetectorMatches) const

1302

{

1303

string locParticleBranchName = "ChargedHypo";

1304

1305

//ASSOCIATED OBJECTS

1306

const DTrackTimeBased* locTrackTimeBased = NULL__null;

1307

locChargedTrackHypothesis->GetSingleT(locTrackTimeBased);

1308

1309

const DBCALShower* locBCALShower = NULL__null;

1310

if(locChargedTrackHypothesis->Get_BCALShowerMatchParams() != NULL__null)

1311

locBCALShower = locChargedTrackHypothesis->Get_BCALShowerMatchParams()->dBCALShower;

1312

1313

const DFCALShower* locFCALShower = NULL__null;

1314

if(locChargedTrackHypothesis->Get_FCALShowerMatchParams() != NULL__null)

1315

locFCALShower = locChargedTrackHypothesis->Get_FCALShowerMatchParams()->dFCALShower;

1316

1317

//IDENTIFIERS

1318

Fill_FundamentalData<Int_t>(locTree, locParticleBranchName, "TrackID", locChargedTrackHypothesis->candidateid, locArrayIndex);

1319

Fill_FundamentalData<Int_t>(locTree, locParticleBranchName, "PID", PDGtype(locChargedTrackHypothesis->PID()), locArrayIndex);

1320

1321

//MATCHING

1322

if(locMCThrownMatching != NULL__null)

1323

{

1324

Int_t locThrownIndex = -1;

1325

double locMatchFOM = 0.0;

1326

const DMCThrown* locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locChargedTrackHypothesis, locMatchFOM);

1327

if(locMCThrown != NULL__null)

1328

locThrownIndex = locThrownIndexMap.find(locMCThrown)->second;

1329

Fill_FundamentalData<Int_t>(locTree, locParticleBranchName, "ThrownIndex", locThrownIndex, locArrayIndex);

1330

}

1331

1332

//KINEMATICS: MEASURED

1333

DVector3 locPosition = locChargedTrackHypothesis->position();

1334

TLorentzVector locTX4_Measured(locPosition.X(), locPosition.Y(), locPosition.Z(), locChargedTrackHypothesis->time());

1335

Fill_ClonesData<TLorentzVector>(locTree, locParticleBranchName, "X4_Measured", locTX4_Measured, locArrayIndex);

1336

1337

DLorentzVector locDP4 = locChargedTrackHypothesis->lorentzMomentum();

1338

TLorentzVector locP4_Measured(locDP4.Px(), locDP4.Py(), locDP4.Pz(), locDP4.E());

1339

Fill_ClonesData<TLorentzVector>(locTree, locParticleBranchName, "P4_Measured", locP4_Measured, locArrayIndex);

1340

1341

//TRACKING INFO

1342

Fill_FundamentalData<UInt_t>(locTree, locParticleBranchName, "NDF_Tracking", locChargedTrackHypothesis->dNDF_Track, locArrayIndex);

1343

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "ChiSq_Tracking", locChargedTrackHypothesis->dChiSq_Track, locArrayIndex);

1344

Fill_FundamentalData<UInt_t>(locTree, locParticleBranchName, "NDF_DCdEdx", locChargedTrackHypothesis->dNDF_DCdEdx, locArrayIndex);

1345

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "ChiSq_DCdEdx", locChargedTrackHypothesis->dChiSq_DCdEdx, locArrayIndex);

1346

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "dEdx_CDC", locTrackTimeBased->ddEdx_CDC, locArrayIndex);

1347

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "dEdx_FDC", locTrackTimeBased->ddEdx_FDC, locArrayIndex);

1348

1349

//HIT ENERGY

1350

double locTOFdEdx = (locChargedTrackHypothesis->Get_TOFHitMatchParams() != NULL__null) ? locChargedTrackHypothesis->Get_TOFHitMatchParams()->dEdx : 0.0;

1351

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "dEdx_TOF", locTOFdEdx, locArrayIndex);

1352

double locSCdEdx = (locChargedTrackHypothesis->Get_SCHitMatchParams() != NULL__null) ? locChargedTrackHypothesis->Get_SCHitMatchParams()->dEdx : 0.0;

1353

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "dEdx_ST", locSCdEdx, locArrayIndex);

1354

double locBCALEnergy = (locBCALShower != NULL__null) ? locBCALShower->E : 0.0;

1355

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "Energy_BCAL", locBCALEnergy, locArrayIndex);

1356

double locFCALEnergy = (locFCALShower != NULL__null) ? locFCALShower->getEnergy() : 0.0;

1357

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "Energy_FCAL", locFCALEnergy, locArrayIndex);

1358

1359

//TIMING INFO

1360

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "HitTime", locChargedTrackHypothesis->t1(), locArrayIndex);

1361

double locStartTimeError = locChargedTrackHypothesis->t0_err();

1362

double locRFDeltaTVariance = (locChargedTrackHypothesis->errorMatrix())(6, 6) + locStartTimeError*locStartTimeError;

1363

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "RFDeltaTVar", locRFDeltaTVariance, locArrayIndex);

1364

1365

//SHOWER MATCHING: BCAL

1366

double locTrackBCAL_DeltaPhi = 999.0, locTrackBCAL_DeltaZ = 999.0;

1367

if(locChargedTrackHypothesis->Get_BCALShowerMatchParams() != NULL__null)

1368

{

1369

locTrackBCAL_DeltaPhi = locChargedTrackHypothesis->Get_BCALShowerMatchParams()->dDeltaPhiToShower;

1370

locTrackBCAL_DeltaZ = locChargedTrackHypothesis->Get_BCALShowerMatchParams()->dDeltaZToShower;

1371

}

1372

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "TrackBCAL_DeltaPhi", locTrackBCAL_DeltaPhi, locArrayIndex);

1373

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "TrackBCAL_DeltaZ", locTrackBCAL_DeltaZ, locArrayIndex);

1374

1375

//SHOWER MATCHING: FCAL

1376

double locDOCAToShower_FCAL = 999.0;

1377

if(locChargedTrackHypothesis->Get_FCALShowerMatchParams() != NULL__null)

1378

locDOCAToShower_FCAL = locChargedTrackHypothesis->Get_FCALShowerMatchParams()->dDOCAToShower;

1379

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "TrackFCAL_DOCA", locDOCAToShower_FCAL, locArrayIndex);

1380

}

1381

1382

void DEventWriterROOT::Fill_NeutralShower(TTree* locTree, unsigned int locArrayIndex, const DNeutralParticleHypothesis* locPhotonHypothesis, const DMCThrownMatching* locMCThrownMatching, const map<const DMCThrown*, unsigned int>& locThrownIndexMap, const DDetectorMatches* locDetectorMatches) const

1383

{

1384

string locParticleBranchName = "NeutralShower";

1385

const DNeutralShower* locNeutralShower = NULL__null;

1386

locPhotonHypothesis->GetSingle(locNeutralShower);

1387

1388

//ASSOCIATED OBJECTS

1389

const DBCALShower* locBCALShower = NULL__null;

1390

locNeutralShower->GetSingle(locBCALShower);

1391

const DFCALShower* locFCALShower = NULL__null;

1392

locNeutralShower->GetSingle(locFCALShower);

1393

1394

//MATCHING

1395

if(locMCThrownMatching != NULL__null)

1396

{

1397

Int_t locThrownIndex = -1;

1398

double locMatchFOM = 0.0;

1399

const DMCThrown* locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locPhotonHypothesis, locMatchFOM);

1400

if(locMCThrown != NULL__null)

1401

locThrownIndex = locThrownIndexMap.find(locMCThrown)->second;

1402

Fill_FundamentalData<Int_t>(locTree, locParticleBranchName, "ThrownIndex", locThrownIndex, locArrayIndex);

1403

}

1404

1405

//SHOWER ENERGY

1406

double locBCALEnergy = (locNeutralShower->dDetectorSystem == SYS_BCAL) ? locNeutralShower->dEnergy : 0.0;

1407

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "Energy_BCAL", locBCALEnergy, locArrayIndex);

1408

double locFCALEnergy = (locNeutralShower->dDetectorSystem == SYS_FCAL) ? locNeutralShower->dEnergy : 0.0;

1409

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "Energy_FCAL", locFCALEnergy, locArrayIndex);

1410

1411

//SHOWER POSITION

1412

DLorentzVector locHitDX4 = locNeutralShower->dSpacetimeVertex;

1413

TLorentzVector locTX4_Shower(locHitDX4.X(), locHitDX4.Y(), locHitDX4.Z(), locHitDX4.T());

1414

Fill_ClonesData<TLorentzVector>(locTree, locParticleBranchName, "X4_Shower", locTX4_Shower, locArrayIndex);

1415

1416

//Track DOCA to Shower - BCAL

1417

double locNearestTrackBCALDeltaPhi = 999.0, locNearestTrackBCALDeltaZ = 999.0;

1418

if(locBCALShower != NULL__null)

1419

{

1420

if(!locDetectorMatches->Get_DistanceToNearestTrack(locBCALShower, locNearestTrackBCALDeltaPhi, locNearestTrackBCALDeltaZ))

1421

{

1422

locNearestTrackBCALDeltaPhi = 999.0;

1423

locNearestTrackBCALDeltaZ = 999.0;

1424

}

1425

else if((locNearestTrackBCALDeltaPhi > 999.0) || (locNearestTrackBCALDeltaZ > 999.0))

1426

{

1427

locNearestTrackBCALDeltaPhi = 999.0;

1428

locNearestTrackBCALDeltaZ = 999.0;

1429

}

1430

}

1431

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "TrackBCAL_DeltaPhi", locNearestTrackBCALDeltaPhi, locArrayIndex);

1432

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "TrackBCAL_DeltaZ", locNearestTrackBCALDeltaZ, locArrayIndex);

1433

1434

//Track DOCA to Shower - FCAL

1435

double locDistanceToNearestTrack_FCAL = 999.0;

1436

if(locFCALShower != NULL__null)

1437

{

1438

if(!locDetectorMatches->Get_DistanceToNearestTrack(locFCALShower, locDistanceToNearestTrack_FCAL))

1439

locDistanceToNearestTrack_FCAL = 999.0;

1440

if(locDistanceToNearestTrack_FCAL > 999.0)

1441

locDistanceToNearestTrack_FCAL = 999.0;

1442

}

1443

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "TrackFCAL_DOCA", locDistanceToNearestTrack_FCAL, locArrayIndex);

1444

1445

//PHOTON PID INFO

1446

double locStartTimeError = locPhotonHypothesis->t0_err();

1447

double PhotonRFDeltaTVar = (locPhotonHypothesis->errorMatrix())(6, 6) + locStartTimeError*locStartTimeError;

1448

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "PhotonRFDeltaTVar", PhotonRFDeltaTVar, locArrayIndex);

1449

}

1450

1451

void DEventWriterROOT::Fill_ComboData(TTree* locTree, const DParticleCombo* locParticleCombo, unsigned int locComboIndex, const map<string, map<oid_t, int> >& locObjectToArrayIndexMap) const

1452

{

1453

//MAIN CLASSES

1454

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

1455

const DKinFitResults* locKinFitResults = locParticleCombo->Get_KinFitResults();

1456

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

1457

1458

//IS COMBO CUT

1459

Fill_FundamentalData<Bool_t>(locTree, "IsComboCut", kFALSE, locComboIndex);

1460

1461

//RF INFO

1462

double locRFTime = (locEventRFBunch != NULL__null) ? locEventRFBunch->dTime : numeric_limits<double>::quiet_NaN();

1463

Fill_FundamentalData<Float_t>(locTree, "RFTime_Measured", locRFTime, locComboIndex);

1464

1465

//KINFIT INFO

1466

DKinFitType locKinFitType = locReaction->Get_KinFitType();

1467

bool locKinFitFlag = (locKinFitType != d_NoFit);

1468

if(locKinFitFlag)

1469

{

1470

if(locKinFitResults != NULL__null)

1471

{

1472

Fill_FundamentalData<Float_t>(locTree, "ChiSq_KinFit", locKinFitResults->Get_ChiSq(), locComboIndex);

1473

Fill_FundamentalData<UInt_t>(locTree, "NDF_KinFit", locKinFitResults->Get_NDF(), locComboIndex);

1474

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

1475

{

1476

double locRFTime_KinFit = -9.9E9; //NOT IMPLEMENTED YET

1477

Fill_FundamentalData<Float_t>(locTree, "RFTime_KinFit", locRFTime_KinFit, locComboIndex);

1478

}

1479

}

1480

else

1481

{

1482

Fill_FundamentalData<Float_t>(locTree, "ChiSq_KinFit", 0.0, locComboIndex);

1483

Fill_FundamentalData<UInt_t>(locTree, "NDF_KinFit", 0, locComboIndex);

1484

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

1485

Fill_FundamentalData<Float_t>(locTree, "RFTime_KinFit", -9.9E9, locComboIndex);

1486

}

1487

}

1488

1489

//STEP DATA

1490

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

1491

Fill_ComboStepData(locTree, locParticleCombo, loc_i, locComboIndex, locKinFitType, locObjectToArrayIndexMap);

1492

}

1493

1494

void DEventWriterROOT::Fill_ComboStepData(TTree* locTree, const DParticleCombo* locParticleCombo, unsigned int locStepIndex, unsigned int locComboIndex, DKinFitType locKinFitType, const map<string, map<oid_t, int> >& locObjectToArrayIndexMap) const

1495

{

1496

TList* locUserInfo = locTree->GetUserInfo();

1497

TMap* locPositionToNameMap = (TMap*)locUserInfo->FindObject("PositionToNameMap");

1498

1499

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

1500

DLorentzVector locStepX4 = locParticleComboStep->Get_SpacetimeVertex();

1501

TLorentzVector locStepTX4(locStepX4.X(), locStepX4.Y(), locStepX4.Z(), locStepX4.T());

1502

1503

//beam & production vertex

1504

Particle_t locInitialPID = locParticleComboStep->Get_InitialParticleID();

1505

const DKinematicData* locInitialParticle = locParticleComboStep->Get_InitialParticle();

1506

const DBeamPhoton* locBeamPhoton = dynamic_cast<const DBeamPhoton*>(locInitialParticle);

1507

if(locBeamPhoton != NULL__null)

1508

{

1509

const DKinematicData* locInitParticleMeasured = locParticleComboStep->Get_InitialParticle_Measured();

1510

const DBeamPhoton* locMeasuredBeamPhoton = dynamic_cast<const DBeamPhoton*>(locInitParticleMeasured);

1511

1512

int locBeamIndex = locObjectToArrayIndexMap.find("DBeamPhoton")->second.find(locMeasuredBeamPhoton->id)->second;

1513

Fill_ComboBeamData(locTree, locComboIndex, locMeasuredBeamPhoton, locBeamPhoton, locBeamIndex, locKinFitType);

1514

}

1515

else //decaying

1516

{

1517

string locParticleBranchName = string("Decaying") + Convert_ToBranchName(ParticleType(locInitialPID));

1518

if(IsDetachedVertex(locInitialPID))

1519

Fill_ClonesData<TLorentzVector>(locTree, locParticleBranchName, "X4", locStepTX4, locComboIndex);

1520

if(IsFixedMass(locInitialPID) && ((locKinFitType == d_P4Fit) || (locKinFitType == d_P4AndVertexFit) || (locKinFitType == d_P4AndSpacetimeFit)))

1521

{

1522

TLorentzVector locDecayP4;

1523

if(locInitialParticle == NULL__null)

1524

{

1525

//fit failed to converge, calc from other particles

1526

DLorentzVector locDecayDP4 = dAnalysisUtilities->Calc_FinalStateP4(locParticleCombo, locStepIndex, false);

1527

locDecayDP4.SetE(sqrt(locDecayDP4.Vect().Mag2() + ParticleMass(locInitialPID)*ParticleMass(locInitialPID)));

1528

locDecayP4.SetPxPyPzE(locDecayDP4.Px(), locDecayDP4.Py(), locDecayDP4.Pz(), locDecayDP4.E());

1529

}

1530

else

1531

locDecayP4.SetPxPyPzE(locInitialParticle->momentum().X(), locInitialParticle->momentum().Y(), locInitialParticle->momentum().Z(), locInitialParticle->energy());

1532

Fill_ClonesData<TLorentzVector>(locTree, locParticleBranchName, "P4_KinFit", locDecayP4, locComboIndex);

1533

}

1534

}

1535

1536

//final state particles

1537

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

1538

{

1539

Particle_t locPID = locParticleComboStep->Get_FinalParticleID(loc_i);

1540

const DKinematicData* locKinematicData = locParticleComboStep->Get_FinalParticle(loc_i);

1541

const DKinematicData* locKinematicData_Measured = locParticleComboStep->Get_FinalParticle_Measured(loc_i);

1542

1543

//decaying particle

1544

if(locParticleComboStep->Is_FinalParticleDecaying(loc_i))

1545

continue;

1546

1547

//missing particle

1548

if(locParticleComboStep->Is_FinalParticleMissing(loc_i))

1549

{

1550

string locParticleBranchName = string("Missing") + Convert_ToBranchName(ParticleType(locPID));

1551

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

1552

{

1553

TLorentzVector locMissingP4;

1554

if(locKinematicData == NULL__null)

1555

{

1556

//fit failed to converge, calc from other particles

1557

DLorentzVector locMissingDP4 = dAnalysisUtilities->Calc_MissingP4(locParticleCombo, false);

1558

locMissingDP4.SetE(sqrt(locMissingDP4.Vect().Mag2() + ParticleMass(locPID)*ParticleMass(locPID)));

1559

locMissingP4.SetPxPyPzE(locMissingDP4.Px(), locMissingDP4.Py(), locMissingDP4.Pz(), locMissingDP4.E());

1560

}

1561

else

1562

locMissingP4.SetPxPyPzE(locKinematicData->momentum().X(), locKinematicData->momentum().Y(), locKinematicData->momentum().Z(), locKinematicData->energy());

1563

1564

Fill_ClonesData<TLorentzVector>(locTree, locParticleBranchName, "P4_KinFit", locMissingP4, locComboIndex);

1565

}

1566

continue;

1567

}

1568

1569

//get the branch name

1570

ostringstream locPositionStream;

1571

locPositionStream << locStepIndex << "_" << loc_i;

1572

TObjString* locObjString = (TObjString*)locPositionToNameMap->GetValue(locPositionStream.str().c_str());

1573

string locParticleBranchName = (const char*)(locObjString->GetString());

1574

1575

//fill the data

1576

if(ParticleCharge(locPID) == 0)

1577

{

1578

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

1579

const DNeutralParticleHypothesis* locMeasuredNeutralHypo = dynamic_cast<const DNeutralParticleHypothesis*>(locKinematicData_Measured);

1580

const DNeutralShower* locNeutralShower = dynamic_cast<const DNeutralShower*>(locParticleComboStep->Get_FinalParticle_SourceObject(loc_i));

1581

int locShowerIndex = locObjectToArrayIndexMap.find("DNeutralShower")->second.find(locNeutralShower->dShowerID)->second;

1582

Fill_ComboNeutralData(locTree, locComboIndex, locParticleBranchName, locMeasuredNeutralHypo, locNeutralHypo, locShowerIndex, locKinFitType);

1583

}

1584

else

1585

{

1586

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

1587

const DChargedTrackHypothesis* locMeasuredChargedHypo = dynamic_cast<const DChargedTrackHypothesis*>(locKinematicData_Measured);

1588

1589

//find "ChargedIndex"

1590

const map<oid_t, int>& locObjectIDMap = locObjectToArrayIndexMap.find("DChargedTrackHypothesis")->second;

1591

vector<const DChargedTrackHypothesis*> locAssociatedChargedHypos;

1592

locKinematicData_Measured->Get(locAssociatedChargedHypos);

1593

int locChargedIndex = -1;

1594

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

1595

{

1596

map<oid_t, int>::const_iterator locIDMapIterator = locObjectIDMap.find(locAssociatedChargedHypos[loc_j]->id);

1597

if(locIDMapIterator == locObjectIDMap.end())

1598

continue;

1599

locChargedIndex = locIDMapIterator->second;

1600

break;

1601

}

1602

if(locChargedIndex == -1)

1603

locChargedIndex = locObjectIDMap.find(locMeasuredChargedHypo->id)->second;

1604

1605

Fill_ComboChargedData(locTree, locComboIndex, locParticleBranchName, locMeasuredChargedHypo, locChargedHypo, locChargedIndex, locKinFitType);

1606

}

1607

}

1608

}

1609

1610

void DEventWriterROOT::Fill_ComboBeamData(TTree* locTree, unsigned int locComboIndex, const DBeamPhoton* locMeasuredBeamPhoton, const DBeamPhoton* locBeamPhoton, unsigned int locBeamIndex, DKinFitType locKinFitType) const

1611

{

1612

string locParticleBranchName = "ComboBeam";

1613

1614

//IDENTIFIER

1615

Fill_FundamentalData<Int_t>(locTree, locParticleBranchName, "BeamIndex", locBeamIndex, locComboIndex);

1616

1617

//MEASURED INTERACTION SPACETIME

1618

DVector3 locPosition = locMeasuredBeamPhoton->position();

1619

TLorentzVector locX4_Measured(locPosition.X(), locPosition.Y(), locPosition.Z(), locMeasuredBeamPhoton->time());

1620

Fill_ClonesData<TLorentzVector>(locTree, locParticleBranchName, "X4_Measured", locX4_Measured, locComboIndex);

1621

1622

//KINEMATICS: KINFIT

1623

if(locKinFitType != d_NoFit)

1624

{

1625

if(locKinFitType != d_P4Fit)

1626

{

1627

DVector3 locPosition = locBeamPhoton->position();

1628

TLorentzVector locX4_KinFit(locPosition.X(), locPosition.Y(), locPosition.Z(), locBeamPhoton->time());

1629

Fill_ClonesData<TLorentzVector>(locTree, locParticleBranchName, "X4_KinFit", locX4_KinFit, locComboIndex);

1630

}

1631

1632

//if charged, bends in b-field, update p4 when vertex changes

1633

if(((locKinFitType != d_VertexFit) && (locKinFitType != d_SpacetimeFit)) || (ParticleCharge(locMeasuredBeamPhoton->PID()) != 0))

1634

{

1635

DLorentzVector locDP4 = locBeamPhoton->lorentzMomentum();

1636

TLorentzVector locP4_KinFit(locDP4.Px(), locDP4.Py(), locDP4.Pz(), locDP4.E());

1637

Fill_ClonesData<TLorentzVector>(locTree, locParticleBranchName, "P4_KinFit", locP4_KinFit, locComboIndex);

1638

}

1639

}

1640

}

1641

1642

void DEventWriterROOT::Fill_ComboChargedData(TTree* locTree, unsigned int locComboIndex, string locParticleBranchName, const DChargedTrackHypothesis* locMeasuredChargedHypo, const DChargedTrackHypothesis* locChargedHypo, unsigned int locChargedIndex, DKinFitType locKinFitType) const

1643

{

1644

//IDENTIFIER

1645

Fill_FundamentalData<Int_t>(locTree, locParticleBranchName, "ChargedIndex", locChargedIndex, locComboIndex);

1646

1647

//MEASURED PID INFO

1648

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "Beta_Timing_Measured", locMeasuredChargedHypo->measuredBeta(), locComboIndex);

1649

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "ChiSq_Timing_Measured", locMeasuredChargedHypo->dChiSq_Timing, locComboIndex);

1650

Fill_FundamentalData<UInt_t>(locTree, locParticleBranchName, "NDF_Timing", locMeasuredChargedHypo->dNDF_Timing, locComboIndex);

1651

1652

//KINFIT

1653

if(locKinFitType != d_NoFit)

1654

{

1655

//KINEMATICS

1656

if(locKinFitType != d_P4Fit)

1657

{

1658

DVector3 locPosition = locChargedHypo->position();

1659

TLorentzVector locX4_KinFit(locPosition.X(), locPosition.Y(), locPosition.Z(), locChargedHypo->time());

1660

Fill_ClonesData<TLorentzVector>(locTree, locParticleBranchName, "X4_KinFit", locX4_KinFit, locComboIndex);

1661

}

1662

1663

//update even if vertex-only fit, because charged momentum propagated through b-field

1664

DLorentzVector locDP4 = locChargedHypo->lorentzMomentum();

1665

TLorentzVector locP4_KinFit(locDP4.Px(), locDP4.Py(), locDP4.Pz(), locDP4.E());

1666

Fill_ClonesData<TLorentzVector>(locTree, locParticleBranchName, "P4_KinFit", locP4_KinFit, locComboIndex);

1667

1668

//PID INFO

1669

if(locKinFitType != d_P4Fit)

1670

{

1671

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "Beta_Timing_KinFit", locChargedHypo->measuredBeta(), locComboIndex);

1672

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "ChiSq_Timing_KinFit", locChargedHypo->dChiSq_Timing, locComboIndex);

1673

}

1674

}

1675

}

1676

1677

void DEventWriterROOT::Fill_ComboNeutralData(TTree* locTree, unsigned int locComboIndex, string locParticleBranchName, const DNeutralParticleHypothesis* locMeasuredNeutralHypo, const DNeutralParticleHypothesis* locNeutralHypo, unsigned int locShowerIndex, DKinFitType locKinFitType) const

1678

{

1679

//IDENTIFIER

1680

Fill_FundamentalData<Int_t>(locTree, locParticleBranchName, "ShowerIndex", locShowerIndex, locComboIndex);

1681

1682

//KINEMATICS: MEASURED

1683

DVector3 locPosition = locMeasuredNeutralHypo->position();

1684

TLorentzVector locX4_Measured(locPosition.X(), locPosition.Y(), locPosition.Z(), locMeasuredNeutralHypo->time());

1685

Fill_ClonesData<TLorentzVector>(locTree, locParticleBranchName, "X4_Measured", locX4_Measured, locComboIndex);

1686

1687

DLorentzVector locDP4 = locMeasuredNeutralHypo->lorentzMomentum();

1688

TLorentzVector locP4_Measured(locDP4.Px(), locDP4.Py(), locDP4.Pz(), locDP4.E());

1689

Fill_ClonesData<TLorentzVector>(locTree, locParticleBranchName, "P4_Measured", locP4_Measured, locComboIndex);

1690

1691

//MEASURED PID INFO

1692

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "Beta_Timing_Measured", locMeasuredNeutralHypo->measuredBeta(), locComboIndex);

1693

if(locParticleBranchName.substr(0, 6) == "Photon")

1694

{

1695

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "ChiSq_Timing_Measured", locMeasuredNeutralHypo->dChiSq, locComboIndex);

1696

Fill_FundamentalData<UInt_t>(locTree, locParticleBranchName, "NDF_Timing", locMeasuredNeutralHypo->dNDF, locComboIndex);

1697

}

1698

1699

//KINFIT

1700

if(locKinFitType != d_NoFit)

1701

{

1702

//KINEMATICS

1703

if(locKinFitType != d_P4Fit)

1704

{

1705

DVector3 locPosition = locNeutralHypo->position();

1706

TLorentzVector locX4_KinFit(locPosition.X(), locPosition.Y(), locPosition.Z(), locNeutralHypo->time());

1707

Fill_ClonesData<TLorentzVector>(locTree, locParticleBranchName, "X4_KinFit", locX4_KinFit, locComboIndex);

1708

}

1709

1710

//update even if vertex-only fit, because neutral momentum defined by vertex

1711

DLorentzVector locDP4 = locNeutralHypo->lorentzMomentum();

1712

TLorentzVector locP4_KinFit(locDP4.Px(), locDP4.Py(), locDP4.Pz(), locDP4.E());

1713

Fill_ClonesData<TLorentzVector>(locTree, locParticleBranchName, "P4_KinFit", locP4_KinFit, locComboIndex);

1714

1715

//PID INFO

1716

if(locKinFitType != d_P4Fit)

1717

{

1718

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "Beta_Timing_KinFit", locNeutralHypo->measuredBeta(), locComboIndex);

1719

if(locParticleBranchName.substr(0, 6) == "Photon")

1720

Fill_FundamentalData<Float_t>(locTree, locParticleBranchName, "ChiSq_Timing_KinFit", locNeutralHypo->dChiSq, locComboIndex);

1721

}

1722

}

1723

}

1724