| File: | /home/sdobbs/work/clang/halld_recon/src/libraries/ANALYSIS/DCutActions.cc |
| Warning: | line 1302, column 14 Value stored to 'locKinFitType' during its initialization is never read |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
| 1 | #ifdef VTRACE |
| 2 | #include "vt_user.h" |
| 3 | #endif |
| 4 | |
| 5 | #include "ANALYSIS/DCutActions.h" |
| 6 | |
| 7 | void DCutAction_MinTrackHits::Initialize(JEventLoop* locEventLoop) |
| 8 | { |
| 9 | Run_Update(locEventLoop); |
| 10 | } |
| 11 | |
| 12 | void DCutAction_MinTrackHits::Run_Update(JEventLoop* locEventLoop) |
| 13 | { |
| 14 | locEventLoop->GetSingle(dParticleID); |
| 15 | } |
| 16 | |
| 17 | |
| 18 | string DCutAction_MinTrackHits::Get_ActionName(void) const |
| 19 | { |
| 20 | ostringstream locStream; |
| 21 | locStream << DAnalysisAction::Get_ActionName() << "_" << dMinTrackHits; |
| 22 | return locStream.str(); |
| 23 | } |
| 24 | |
| 25 | bool DCutAction_MinTrackHits::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 26 | { |
| 27 | auto locParticles = locParticleCombo->Get_FinalParticles_Measured(Get_Reaction(), d_Charged); |
| 28 | for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) |
| 29 | { |
| 30 | const DChargedTrackHypothesis* locChargedTrackHypothesis = static_cast<const DChargedTrackHypothesis*>(locParticles[loc_i]); |
| 31 | auto locTrackTimeBased = locChargedTrackHypothesis->Get_TrackTimeBased(); |
| 32 | |
| 33 | set<int> locCDCRings, locFDCPlanes; |
| 34 | dParticleID->Get_CDCRings(locTrackTimeBased->dCDCRings, locCDCRings); |
| 35 | dParticleID->Get_FDCPlanes(locTrackTimeBased->dFDCPlanes, locFDCPlanes); |
| 36 | unsigned int locNumTrackHits = locCDCRings.size() + locFDCPlanes.size(); |
| 37 | if(locNumTrackHits < dMinTrackHits) |
| 38 | return false; |
| 39 | } |
| 40 | return true; |
| 41 | } |
| 42 | |
| 43 | string DCutAction_ThrownTopology::Get_ActionName(void) const |
| 44 | { |
| 45 | ostringstream locStream; |
| 46 | locStream << DAnalysisAction::Get_ActionName() << "_" << dExclusiveMatchFlag; |
| 47 | return locStream.str(); |
| 48 | } |
| 49 | |
| 50 | void DCutAction_ThrownTopology::Initialize(JEventLoop* locEventLoop) |
| 51 | { |
| 52 | Run_Update(locEventLoop); |
| 53 | } |
| 54 | |
| 55 | void DCutAction_ThrownTopology::Run_Update(JEventLoop* locEventLoop) |
| 56 | { |
| 57 | locEventLoop->GetSingle(dAnalysisUtilities); |
| 58 | } |
| 59 | |
| 60 | bool DCutAction_ThrownTopology::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 61 | { |
| 62 | return dAnalysisUtilities->Check_ThrownsMatchReaction(locEventLoop, Get_Reaction(), dExclusiveMatchFlag); |
| 63 | } |
| 64 | |
| 65 | bool DCutAction_AllTracksHaveDetectorMatch::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 66 | { |
| 67 | auto locParticles = locParticleCombo->Get_FinalParticles_Measured(Get_Reaction(), d_Charged); |
| 68 | for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) |
| 69 | { |
| 70 | const DChargedTrackHypothesis* locChargedTrackHypothesis = static_cast<const DChargedTrackHypothesis*>(locParticles[loc_i]); |
| 71 | if(locChargedTrackHypothesis->Get_SCHitMatchParams() != NULL__null) |
| 72 | continue; |
| 73 | if(locChargedTrackHypothesis->Get_TOFHitMatchParams() != NULL__null) |
| 74 | continue; |
| 75 | if(locChargedTrackHypothesis->Get_BCALShowerMatchParams() != NULL__null) |
| 76 | continue; |
| 77 | if(locChargedTrackHypothesis->Get_FCALShowerMatchParams() != NULL__null) |
| 78 | continue; |
| 79 | return false; |
| 80 | } |
| 81 | return true; |
| 82 | } |
| 83 | |
| 84 | string DCutAction_PIDFOM::Get_ActionName(void) const |
| 85 | { |
| 86 | ostringstream locStream; |
| 87 | locStream << DAnalysisAction::Get_ActionName() << "_" << dMinimumConfidenceLevel; |
| 88 | return locStream.str(); |
| 89 | } |
| 90 | |
| 91 | bool DCutAction_PIDFOM::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 92 | { |
| 93 | auto locSteps = locParticleCombo->Get_ParticleComboSteps(); |
| 94 | for(size_t loc_i = 0; loc_i < locSteps.size(); ++loc_i) |
| 95 | { |
| 96 | if((dStepPID != Unknown) && (Get_Reaction()->Get_ReactionStep(loc_i)->Get_InitialPID() != dStepPID)) |
| 97 | continue; |
| 98 | auto locParticles = locSteps[loc_i]->Get_FinalParticles_Measured(Get_Reaction()->Get_ReactionStep(loc_i), d_AllCharges); |
| 99 | for(size_t loc_j = 0; loc_j < locParticles.size(); ++loc_j) |
| 100 | { |
| 101 | if((locParticles[loc_j]->PID() != dParticleID) && (dParticleID != Unknown)) |
| 102 | continue; |
| 103 | if(ParticleCharge(dParticleID) == 0) |
| 104 | { |
| 105 | const DNeutralParticleHypothesis* locNeutralParticleHypothesis = static_cast<const DNeutralParticleHypothesis*>(locParticles[loc_i]); |
| 106 | if((locNeutralParticleHypothesis->Get_FOM() < dMinimumConfidenceLevel) && (locNeutralParticleHypothesis->Get_NDF() > 0)) |
| 107 | return false; |
| 108 | if(dCutNDFZeroFlag && (locNeutralParticleHypothesis->Get_NDF() == 0)) |
| 109 | return false; |
| 110 | } |
| 111 | else |
| 112 | { |
| 113 | const DChargedTrackHypothesis* locChargedTrackHypothesis = static_cast<const DChargedTrackHypothesis*>(locParticles[loc_i]); |
| 114 | if((locChargedTrackHypothesis->Get_FOM() < dMinimumConfidenceLevel) && (locChargedTrackHypothesis->Get_NDF() > 0)) |
| 115 | return false; |
| 116 | if(dCutNDFZeroFlag && (locChargedTrackHypothesis->Get_NDF() == 0)) |
| 117 | return false; |
| 118 | } |
| 119 | } |
| 120 | } |
| 121 | return true; |
| 122 | } |
| 123 | |
| 124 | string DCutAction_EachPIDFOM::Get_ActionName(void) const |
| 125 | { |
| 126 | ostringstream locStream; |
| 127 | locStream << DAnalysisAction::Get_ActionName() << "_" << dMinimumConfidenceLevel; |
| 128 | return locStream.str(); |
| 129 | } |
| 130 | |
| 131 | bool DCutAction_EachPIDFOM::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 132 | { |
| 133 | auto locParticles = locParticleCombo->Get_FinalParticles_Measured(Get_Reaction(), d_AllCharges); |
| 134 | for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) |
| 135 | { |
| 136 | if(ParticleCharge(locParticles[loc_i]->PID()) == 0) |
| 137 | { |
| 138 | const DNeutralParticleHypothesis* locNeutralParticleHypothesis = static_cast<const DNeutralParticleHypothesis*>(locParticles[loc_i]); |
| 139 | if(dCutNDFZeroFlag && (locNeutralParticleHypothesis->Get_NDF() == 0)) |
| 140 | return false; |
| 141 | if((locNeutralParticleHypothesis->Get_FOM() < dMinimumConfidenceLevel) && (locNeutralParticleHypothesis->Get_NDF() > 0)) |
| 142 | return false; |
| 143 | } |
| 144 | else |
| 145 | { |
| 146 | const DChargedTrackHypothesis* locChargedTrackHypothesis = static_cast<const DChargedTrackHypothesis*>(locParticles[loc_i]); |
| 147 | if(dCutNDFZeroFlag && (locChargedTrackHypothesis->Get_NDF() == 0)) |
| 148 | return false; |
| 149 | if((locChargedTrackHypothesis->Get_FOM() < dMinimumConfidenceLevel) && (locChargedTrackHypothesis->Get_NDF() > 0)) |
| 150 | return false; |
| 151 | } |
| 152 | } |
| 153 | return true; |
| 154 | } |
| 155 | |
| 156 | string DCutAction_CombinedPIDFOM::Get_ActionName(void) const |
| 157 | { |
| 158 | ostringstream locStream; |
| 159 | locStream << DAnalysisAction::Get_ActionName() << "_" << dMinimumConfidenceLevel; |
| 160 | return locStream.str(); |
| 161 | } |
| 162 | |
| 163 | bool DCutAction_CombinedPIDFOM::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 164 | { |
| 165 | auto locParticles = locParticleCombo->Get_FinalParticles_Measured(Get_Reaction(), d_AllCharges); |
| 166 | |
| 167 | unsigned int locTotalNDF = 0; |
| 168 | double locTotalChiSq = 0.0; |
| 169 | for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) |
| 170 | { |
| 171 | if(ParticleCharge(locParticles[loc_i]->PID()) == 0) |
| 172 | { |
| 173 | const DNeutralParticleHypothesis* locNeutralParticleHypothesis = static_cast<const DNeutralParticleHypothesis*>(locParticles[loc_i]); |
| 174 | if(dCutNDFZeroFlag && (locNeutralParticleHypothesis->Get_NDF() == 0)) |
| 175 | return false; |
| 176 | locTotalNDF += locNeutralParticleHypothesis->Get_NDF(); |
| 177 | locTotalChiSq += locNeutralParticleHypothesis->Get_ChiSq(); |
| 178 | } |
| 179 | else |
| 180 | { |
| 181 | const DChargedTrackHypothesis* locChargedTrackHypothesis = static_cast<const DChargedTrackHypothesis*>(locParticles[loc_i]); |
| 182 | if(dCutNDFZeroFlag && (locChargedTrackHypothesis->Get_NDF() == 0)) |
| 183 | return false; |
| 184 | locTotalNDF += locChargedTrackHypothesis->Get_NDF(); |
| 185 | locTotalChiSq += locChargedTrackHypothesis->Get_ChiSq(); |
| 186 | } |
| 187 | } |
| 188 | return ((locTotalNDF == 0) ? true : (TMath::Prob(locTotalChiSq, locTotalNDF) >= dMinimumConfidenceLevel)); |
| 189 | } |
| 190 | |
| 191 | string DCutAction_CombinedTrackingFOM::Get_ActionName(void) const |
| 192 | { |
| 193 | ostringstream locStream; |
| 194 | locStream << DAnalysisAction::Get_ActionName() << "_" << dMinimumConfidenceLevel; |
| 195 | return locStream.str(); |
| 196 | } |
| 197 | |
| 198 | bool DCutAction_CombinedTrackingFOM::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 199 | { |
| 200 | unsigned int locTotalNDF = 0; |
| 201 | double locTotalChiSq = 0.0; |
| 202 | |
| 203 | auto locParticles = locParticleCombo->Get_FinalParticles_Measured(Get_Reaction(), d_Charged); |
| 204 | for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) |
| 205 | { |
| 206 | auto locTrackTimeBased = (dynamic_cast<const DChargedTrackHypothesis*>(locParticles[loc_i]))->Get_TrackTimeBased(); |
| 207 | locTotalNDF += locTrackTimeBased->Ndof; |
| 208 | locTotalChiSq += locTrackTimeBased->chisq; |
| 209 | } |
| 210 | |
| 211 | return ((locTotalNDF == 0) ? true : (TMath::Prob(locTotalChiSq, locTotalNDF) >= dMinimumConfidenceLevel)); |
| 212 | } |
| 213 | |
| 214 | string DCutAction_MissingMass::Get_ActionName(void) const |
| 215 | { |
| 216 | ostringstream locStream; |
| 217 | locStream << DAnalysisAction::Get_ActionName() << "_" << dMinimumMissingMass << "_" << dMaximumMissingMass; |
| 218 | return locStream.str(); |
| 219 | } |
| 220 | |
| 221 | void DCutAction_MissingMass::Initialize(JEventLoop* locEventLoop) |
| 222 | { |
| 223 | Run_Update(locEventLoop); |
| 224 | } |
| 225 | |
| 226 | void DCutAction_MissingMass::Run_Update(JEventLoop* locEventLoop) |
| 227 | { |
| 228 | locEventLoop->GetSingle(dAnalysisUtilities); |
| 229 | } |
| 230 | |
| 231 | bool DCutAction_MissingMass::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 232 | { |
| 233 | //build all possible combinations of the included pids |
| 234 | set<set<size_t> > locIndexCombos = dAnalysisUtilities->Build_IndexCombos(Get_Reaction()->Get_ReactionStep(dMissingMassOffOfStepIndex), dMissingMassOffOfPIDs); |
| 235 | |
| 236 | //loop over them: Must fail ALL to fail. if any succeed, return true |
| 237 | set<set<size_t> >::iterator locComboIterator = locIndexCombos.begin(); |
| 238 | for(; locComboIterator != locIndexCombos.end(); ++locComboIterator) |
| 239 | { |
| 240 | DLorentzVector locMissingP4 = dAnalysisUtilities->Calc_MissingP4(Get_Reaction(), locParticleCombo, 0, dMissingMassOffOfStepIndex, *locComboIterator, Get_UseKinFitResultsFlag()); |
| 241 | double locMissingMass = locMissingP4.M(); |
| 242 | if((locMissingMass >= dMinimumMissingMass) && (locMissingMass <= dMaximumMissingMass)) |
| 243 | return true; |
| 244 | } |
| 245 | |
| 246 | return false; //all failed |
| 247 | } |
| 248 | |
| 249 | string DCutAction_MissingMassSquared::Get_ActionName(void) const |
| 250 | { |
| 251 | ostringstream locStream; |
| 252 | locStream << DAnalysisAction::Get_ActionName() << "_" << dMinimumMissingMassSq << "_" << dMaximumMissingMassSq; |
| 253 | return locStream.str(); |
| 254 | } |
| 255 | |
| 256 | void DCutAction_MissingMassSquared::Initialize(JEventLoop* locEventLoop) |
| 257 | { |
| 258 | locEventLoop->GetSingle(dAnalysisUtilities); |
| 259 | } |
| 260 | |
| 261 | void DCutAction_MissingMassSquared::Run_Update(JEventLoop* locEventLoop) |
| 262 | { |
| 263 | locEventLoop->GetSingle(dAnalysisUtilities); |
| 264 | } |
| 265 | |
| 266 | |
| 267 | bool DCutAction_MissingMassSquared::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 268 | { |
| 269 | //build all possible combinations of the included pids |
| 270 | set<set<size_t> > locIndexCombos = dAnalysisUtilities->Build_IndexCombos(Get_Reaction()->Get_ReactionStep(dMissingMassOffOfStepIndex), dMissingMassOffOfPIDs); |
| 271 | |
| 272 | //loop over them: Must fail ALL to fail. if any succeed, return true |
| 273 | set<set<size_t> >::iterator locComboIterator = locIndexCombos.begin(); |
| 274 | for(; locComboIterator != locIndexCombos.end(); ++locComboIterator) |
| 275 | { |
| 276 | DLorentzVector locMissingP4 = dAnalysisUtilities->Calc_MissingP4(Get_Reaction(), locParticleCombo, 0, dMissingMassOffOfStepIndex, *locComboIterator, Get_UseKinFitResultsFlag()); |
| 277 | double locMissingMassSq = locMissingP4.M2(); |
| 278 | if((locMissingMassSq >= dMinimumMissingMassSq) && (locMissingMassSq <= dMaximumMissingMassSq)) |
| 279 | return true; |
| 280 | } |
| 281 | |
| 282 | return false; //all failed |
| 283 | } |
| 284 | |
| 285 | string DCutAction_InvariantMass::Get_ActionName(void) const |
| 286 | { |
| 287 | ostringstream locStream; |
| 288 | locStream << DAnalysisAction::Get_ActionName() << "_" << dInitialPID << "_" << dMinMass << "_" << dMaxMass; |
| 289 | return locStream.str(); |
| 290 | } |
| 291 | |
| 292 | void DCutAction_InvariantMass::Initialize(JEventLoop* locEventLoop) |
| 293 | { |
| 294 | Run_Update(locEventLoop); |
| 295 | } |
| 296 | |
| 297 | void DCutAction_InvariantMass::Run_Update(JEventLoop* locEventLoop) |
| 298 | { |
| 299 | locEventLoop->GetSingle(dAnalysisUtilities); |
| 300 | } |
| 301 | |
| 302 | bool DCutAction_InvariantMass::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 303 | { |
| 304 | for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i) |
| 305 | { |
| 306 | const DReactionStep* locReactionStep = Get_Reaction()->Get_ReactionStep(loc_i); |
| 307 | if((dInitialPID != Unknown) && (locReactionStep->Get_InitialPID() != dInitialPID)) |
| 308 | continue; |
| 309 | if((dStepIndex != -1) && (int(loc_i) != dStepIndex)) |
| 310 | continue; |
| 311 | |
| 312 | //build all possible combinations of the included pids |
| 313 | set<set<size_t> > locIndexCombos = dAnalysisUtilities->Build_IndexCombos(locReactionStep, dToIncludePIDs); |
| 314 | |
| 315 | //loop over them: Must fail ALL to fail. if any succeed, go to the next step |
| 316 | set<set<size_t> >::iterator locComboIterator = locIndexCombos.begin(); |
| 317 | bool locAnyOKFlag = false; |
| 318 | for(; locComboIterator != locIndexCombos.end(); ++locComboIterator) |
| 319 | { |
| 320 | DLorentzVector locFinalStateP4 = dAnalysisUtilities->Calc_FinalStateP4(Get_Reaction(), locParticleCombo, loc_i, *locComboIterator, Get_UseKinFitResultsFlag()); |
| 321 | double locInvariantMass = locFinalStateP4.M(); |
| 322 | if((locInvariantMass > dMaxMass) || (locInvariantMass < dMinMass)) |
| 323 | continue; |
| 324 | locAnyOKFlag = true; |
| 325 | break; |
| 326 | } |
| 327 | if(!locAnyOKFlag) |
| 328 | return false; |
| 329 | } |
| 330 | |
| 331 | return true; |
| 332 | } |
| 333 | |
| 334 | string DCutAction_AllVertexZ::Get_ActionName(void) const |
| 335 | { |
| 336 | ostringstream locStream; |
| 337 | locStream << DAnalysisAction::Get_ActionName() << "_" << dMinVertexZ << "_" << dMaxVertexZ; |
| 338 | return locStream.str(); |
| 339 | } |
| 340 | |
| 341 | bool DCutAction_AllVertexZ::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 342 | { |
| 343 | auto locParticles = locParticleCombo->Get_FinalParticles_Measured(Get_Reaction(), d_Charged); |
| 344 | double locVertexZ; |
| 345 | for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) |
| 346 | { |
| 347 | locVertexZ = locParticles[loc_i]->position().Z(); |
| 348 | if((locVertexZ < dMinVertexZ) || (locVertexZ > dMaxVertexZ)) |
| 349 | return false; |
| 350 | } |
| 351 | return true; |
| 352 | } |
| 353 | |
| 354 | string DCutAction_ProductionVertexZ::Get_ActionName(void) const |
| 355 | { |
| 356 | ostringstream locStream; |
| 357 | locStream << DAnalysisAction::Get_ActionName() << "_" << dMinVertexZ << "_" << dMaxVertexZ; |
| 358 | return locStream.str(); |
| 359 | } |
| 360 | |
| 361 | bool DCutAction_ProductionVertexZ::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 362 | { |
| 363 | const DParticleComboStep* locStep = locParticleCombo->Get_ParticleComboStep(0); |
| 364 | double locVertexZ = locStep->Get_Position().Z(); |
| 365 | return ((locVertexZ >= dMinVertexZ) && (locVertexZ <= dMaxVertexZ)); |
| 366 | } |
| 367 | |
| 368 | string DCutAction_MaxTrackDOCA::Get_ActionName(void) const |
| 369 | { |
| 370 | ostringstream locStream; |
| 371 | locStream << DAnalysisAction::Get_ActionName() << "_" << dMaxTrackDOCA; |
| 372 | return locStream.str(); |
| 373 | } |
| 374 | |
| 375 | void DCutAction_MaxTrackDOCA::Initialize(JEventLoop* locEventLoop) |
| 376 | { |
| 377 | Run_Update(locEventLoop); |
| 378 | } |
| 379 | |
| 380 | void DCutAction_MaxTrackDOCA::Run_Update(JEventLoop* locEventLoop) |
| 381 | { |
| 382 | locEventLoop->GetSingle(dAnalysisUtilities); |
| 383 | } |
| 384 | |
| 385 | bool DCutAction_MaxTrackDOCA::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 386 | { |
| 387 | //should be improved...: the particles at a given vertex may span several steps |
| 388 | auto locSteps = locParticleCombo->Get_ParticleComboSteps(); |
| 389 | for(size_t loc_i = 0; loc_i < locSteps.size(); ++loc_i) |
| 390 | { |
| 391 | if((dInitialPID != Unknown) && (Get_Reaction()->Get_ReactionStep(loc_i)->Get_InitialPID() != dInitialPID)) |
| 392 | continue; |
| 393 | auto locParticles = locSteps[loc_i]->Get_FinalParticles_Measured(Get_Reaction()->Get_ReactionStep(loc_i), d_Charged); |
| 394 | for(size_t loc_j = 0; loc_j < locParticles.size(); ++loc_j) |
| 395 | { |
| 396 | for(size_t loc_k = loc_j + 1; loc_k < locParticles.size(); ++loc_k) |
| 397 | { |
| 398 | auto locDOCA = dAnalysisUtilities->Calc_DOCA(locParticles[loc_j], locParticles[loc_k]); |
| 399 | if(locDOCA > dMaxTrackDOCA) |
| 400 | return false; |
| 401 | } |
| 402 | } |
| 403 | } |
| 404 | return true; |
| 405 | } |
| 406 | |
| 407 | string DCutAction_KinFitFOM::Get_ActionName(void) const |
| 408 | { |
| 409 | ostringstream locStream; |
| 410 | locStream << DAnalysisAction::Get_ActionName() << "_" << dMinimumConfidenceLevel; |
| 411 | return locStream.str(); |
| 412 | } |
| 413 | |
| 414 | bool DCutAction_KinFitFOM::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 415 | { |
| 416 | const DKinFitResults* locKinFitResults = locParticleCombo->Get_KinFitResults(); |
| 417 | if(locKinFitResults == NULL__null) |
| 418 | return false; |
| 419 | return (locKinFitResults->Get_ConfidenceLevel() > dMinimumConfidenceLevel); |
| 420 | } |
| 421 | |
| 422 | string DCutAction_KinFitChiSq::Get_ActionName(void) const |
| 423 | { |
| 424 | ostringstream locStream; |
| 425 | locStream << DAnalysisAction::Get_ActionName() << "_" << dMaximumChiSq; |
| 426 | return locStream.str(); |
| 427 | } |
| 428 | |
| 429 | bool DCutAction_KinFitChiSq::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 430 | { |
| 431 | const DKinFitResults* locKinFitResults = locParticleCombo->Get_KinFitResults(); |
| 432 | if(locKinFitResults == NULL__null) |
| 433 | return false; |
| 434 | return (locKinFitResults->Get_ChiSq()/locKinFitResults->Get_NDF() < dMaximumChiSq); |
| 435 | } |
| 436 | |
| 437 | void DCutAction_BDTSignalCombo::Initialize(JEventLoop* locEventLoop) |
| 438 | { |
| 439 | if(dCutAction_TrueBeamParticle == nullptr) |
| 440 | dCutAction_TrueBeamParticle = new DCutAction_TrueBeamParticle(Get_Reaction()); |
| 441 | dCutAction_TrueBeamParticle->Initialize(locEventLoop); |
| 442 | Run_Update(locEventLoop); |
| 443 | } |
| 444 | |
| 445 | void DCutAction_BDTSignalCombo::Run_Update(JEventLoop* locEventLoop) |
| 446 | { |
| 447 | locEventLoop->GetSingle(dAnalysisUtilities); |
| 448 | } |
| 449 | |
| 450 | bool DCutAction_BDTSignalCombo::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 451 | { |
| 452 | #ifdef VTRACE |
| 453 | VT_TRACER("DCutAction_BDTSignalCombo::Perform_Action()"); |
| 454 | #endif |
| 455 | |
| 456 | vector<const DMCThrownMatching*> locMCThrownMatchingVector; |
| 457 | locEventLoop->Get(locMCThrownMatchingVector); |
| 458 | if(locMCThrownMatchingVector.empty()) |
| 459 | return false; //not a simulated event |
| 460 | const DMCThrownMatching* locMCThrownMatching = locMCThrownMatchingVector[0]; |
| 461 | |
| 462 | //Check DReaction vs thrown (i.e. not combo contents) |
| 463 | if(!dAnalysisUtilities->Check_IsBDTSignalEvent(locEventLoop, Get_Reaction(), dExclusiveMatchFlag, dIncludeDecayingToReactionFlag)) |
| 464 | return false; |
| 465 | |
| 466 | //Do we need to pick the beam photon? If so, look for it |
| 467 | Particle_t locPID = Get_Reaction()->Get_ReactionStep(0)->Get_InitialPID(); |
| 468 | if(locPID == Gamma) |
| 469 | { |
| 470 | if(!(*dCutAction_TrueBeamParticle)(locEventLoop, locParticleCombo)) |
| 471 | return false; //needed the true beam photon, didn't have it |
| 472 | } |
| 473 | |
| 474 | //get & organize throwns |
| 475 | vector<const DMCThrown*> locMCThrowns; |
| 476 | locEventLoop->Get(locMCThrowns); |
| 477 | |
| 478 | map<int, const DMCThrown*> locMCThrownMyIDMap; //map of myid -> thrown |
| 479 | for(size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i) |
| 480 | locMCThrownMyIDMap[locMCThrowns[loc_i]->myid] = locMCThrowns[loc_i]; |
| 481 | |
| 482 | //OK, now need to check and see if the particles have the right PID & the right parent chain |
| 483 | for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i) |
| 484 | { |
| 485 | const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(loc_i); |
| 486 | auto locReactionStep = Get_Reaction()->Get_ReactionStep(loc_i); |
| 487 | |
| 488 | auto locParticles = locParticleComboStep->Get_FinalParticles_Measured(locReactionStep, d_AllCharges); |
| 489 | for(size_t loc_j = 0; loc_j < locParticles.size(); ++loc_j) |
| 490 | { |
| 491 | const DMCThrown* locMCThrown = NULL__null; |
| 492 | double locMatchFOM = 0.0; |
| 493 | if(ParticleCharge(locParticles[loc_j]->PID()) == 0) |
| 494 | { |
| 495 | //check if good neutral & PID |
| 496 | const DNeutralParticleHypothesis* locNeutralParticleHypothesis = static_cast<const DNeutralParticleHypothesis*>(locParticles[loc_j]); |
| 497 | locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locNeutralParticleHypothesis, locMatchFOM); |
| 498 | if((locMCThrown == NULL__null) || (locMatchFOM < dMinThrownMatchFOM)) |
| 499 | return false; //not matched |
| 500 | if(((Particle_t)locMCThrown->type) != locParticles[loc_j]->PID()) |
| 501 | return false; //bad PID |
| 502 | } |
| 503 | else |
| 504 | { |
| 505 | //check if good track & PID |
| 506 | double locMatchFOM = 0.0; |
| 507 | const DChargedTrackHypothesis* locChargedTrackHypothesis = static_cast<const DChargedTrackHypothesis*>(locParticles[loc_j]); |
| 508 | locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locChargedTrackHypothesis, locMatchFOM); |
| 509 | if((locMCThrown == NULL__null) || (locMatchFOM < dMinThrownMatchFOM)) |
| 510 | return false; //not matched |
| 511 | if(((Particle_t)locMCThrown->type) != locParticles[loc_j]->PID()) |
| 512 | return false; //bad PID |
| 513 | } |
| 514 | |
| 515 | //check if parent is correct |
| 516 | auto locParentSearchReactionStep = locReactionStep; |
| 517 | int locParentID = locMCThrown->parentid; |
| 518 | int locSearchStepIndex = loc_i; |
| 519 | |
| 520 | do |
| 521 | { |
| 522 | if(locParentID == -1) |
| 523 | return false; //parent particle is not listed: matched to knock-out particle, is wrong. bail |
| 524 | |
| 525 | if(locParentID == 0) |
| 526 | { |
| 527 | //parent id of 0 is directly produced |
| 528 | Particle_t locInitPID = locParentSearchReactionStep->Get_InitialPID(); |
| 529 | if((locInitPID == phiMeson) || (locInitPID == omega)) |
| 530 | { |
| 531 | //these particles are not constrained: check their parent step instead |
| 532 | locSearchStepIndex = DAnalysis::Get_InitialParticleDecayFromIndices(Get_Reaction(), locSearchStepIndex).first; |
| 533 | locReactionStep = Get_Reaction()->Get_ReactionStep(locSearchStepIndex); |
| 534 | continue; |
| 535 | } |
| 536 | if(locInitPID != Gamma) |
| 537 | return false; //was directly (photo-) produced, but not so in combo: bail |
| 538 | break; //good: this is the only "good" exit point of the do-loop |
| 539 | } |
| 540 | |
| 541 | const DMCThrown* locMCThrownParent = locMCThrownMyIDMap[locParentID]; |
| 542 | Particle_t locPID = locMCThrownParent->PID(); |
| 543 | if((locPID == Unknown) || IsResonance(locPID) || (locPID == omega) || (locPID == phiMeson)) |
| 544 | { |
| 545 | //intermediate (unknown, resonance, phi, or omega) particle: go to its parent |
| 546 | locParentID = locMCThrownParent->parentid; |
| 547 | continue; |
| 548 | } |
| 549 | |
| 550 | if(locPID != locParentSearchReactionStep->Get_InitialPID()) |
| 551 | { |
| 552 | //the true particle was produced from a different parent |
| 553 | if(!dIncludeDecayingToReactionFlag) |
| 554 | return false; //will not consider intermediate decays: bail |
| 555 | |
| 556 | //it could still be BDT signal, if the thrown parent eventually came from the combo parent particle |
| 557 | //treat as an intermediate decaying particle: go to its parent |
| 558 | locParentID = locMCThrownParent->parentid; |
| 559 | continue; |
| 560 | } |
| 561 | |
| 562 | //OK, we've determined that the particle in question decayed from the correct particle. |
| 563 | //HOWEVER, we need to determine whether the PARENT decayed from the correct particle (and on(back)wards until the production step) |
| 564 | locParentID = locMCThrownParent->parentid; |
| 565 | locSearchStepIndex = DAnalysis::Get_InitialParticleDecayFromIndices(Get_Reaction(), locSearchStepIndex).first; |
| 566 | locReactionStep = Get_Reaction()->Get_ReactionStep(locSearchStepIndex); |
| 567 | } |
| 568 | while(true); |
| 569 | } |
| 570 | } |
| 571 | |
| 572 | return true; //we made it! |
| 573 | } |
| 574 | |
| 575 | DCutAction_BDTSignalCombo::~DCutAction_BDTSignalCombo(void) |
| 576 | { |
| 577 | if(dCutAction_TrueBeamParticle != NULL__null) |
| 578 | delete dCutAction_TrueBeamParticle; |
| 579 | } |
| 580 | |
| 581 | void DCutAction_TrueCombo::Initialize(JEventLoop* locEventLoop) |
| 582 | { |
| 583 | if(dCutAction_TrueBeamParticle == nullptr) |
| 584 | dCutAction_TrueBeamParticle = new DCutAction_TrueBeamParticle(Get_Reaction()); |
| 585 | dCutAction_TrueBeamParticle->Initialize(locEventLoop); |
| 586 | if(dCutAction_ThrownTopology == nullptr) |
| 587 | dCutAction_ThrownTopology = new DCutAction_ThrownTopology(Get_Reaction(), dExclusiveMatchFlag); |
| 588 | dCutAction_ThrownTopology->Initialize(locEventLoop); |
| 589 | } |
| 590 | |
| 591 | void DCutAction_TrueCombo::Run_Update(JEventLoop* locEventLoop) |
| 592 | { |
| 593 | dCutAction_TrueBeamParticle->Run_Update(locEventLoop); |
| 594 | dCutAction_ThrownTopology->Run_Update(locEventLoop); |
| 595 | } |
| 596 | |
| 597 | bool DCutAction_TrueCombo::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 598 | { |
| 599 | #ifdef VTRACE |
| 600 | VT_TRACER("DCutAction_TrueCombo::Perform_Action()"); |
| 601 | #endif |
| 602 | |
| 603 | vector<const DMCThrownMatching*> locMCThrownMatchingVector; |
| 604 | locEventLoop->Get(locMCThrownMatchingVector); |
| 605 | if(locMCThrownMatchingVector.empty()) |
| 606 | return false; //not a simulated event |
| 607 | const DMCThrownMatching* locMCThrownMatching = locMCThrownMatchingVector[0]; |
| 608 | |
| 609 | if(!(*dCutAction_ThrownTopology)(locEventLoop, locParticleCombo)) |
| 610 | return false; //not the thrown topology: bail |
| 611 | |
| 612 | //Do we need to pick the beam photon? If so, look for it |
| 613 | if(DAnalysis::Get_IsFirstStepBeam(Get_Reaction())) |
| 614 | { |
| 615 | if(!(*dCutAction_TrueBeamParticle)(locEventLoop, locParticleCombo)) |
| 616 | return false; //needed the true beam photon, didn't have it |
| 617 | } |
| 618 | |
| 619 | //get & organize throwns |
| 620 | vector<const DMCThrown*> locMCThrowns; |
| 621 | locEventLoop->Get(locMCThrowns); |
| 622 | |
| 623 | map<int, const DMCThrown*> locMCThrownMyIDMap; //map of myid -> thrown |
| 624 | for(size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i) |
| 625 | locMCThrownMyIDMap[locMCThrowns[loc_i]->myid] = locMCThrowns[loc_i]; |
| 626 | |
| 627 | //OK, now need to check and see if the particles have the right PID & the right parent chain |
| 628 | for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i) |
| 629 | { |
| 630 | const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(loc_i); |
| 631 | auto locReactionStep = Get_Reaction()->Get_ReactionStep(loc_i); |
| 632 | |
| 633 | auto locParticles = locParticleComboStep->Get_FinalParticles_Measured(locReactionStep, d_AllCharges); |
| 634 | for(size_t loc_j = 0; loc_j < locParticles.size(); ++loc_j) |
| 635 | { |
| 636 | const DMCThrown* locMCThrown = NULL__null; |
| 637 | double locMatchFOM = 0.0; |
| 638 | if(ParticleCharge(locParticles[loc_j]->PID()) == 0) |
| 639 | { |
| 640 | //check if good neutral & PID |
| 641 | const DNeutralParticleHypothesis* locNeutralParticleHypothesis = static_cast<const DNeutralParticleHypothesis*>(locParticles[loc_j]); |
| 642 | locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locNeutralParticleHypothesis, locMatchFOM); |
| 643 | if((locMCThrown == NULL__null) || (locMatchFOM < dMinThrownMatchFOM)) |
| 644 | return false; //not matched |
| 645 | if(((Particle_t)locMCThrown->type) != locParticles[loc_j]->PID()) |
| 646 | return false; //bad PID |
| 647 | } |
| 648 | else |
| 649 | { |
| 650 | //check if good track & PID |
| 651 | double locMatchFOM = 0.0; |
| 652 | const DChargedTrackHypothesis* locChargedTrackHypothesis = static_cast<const DChargedTrackHypothesis*>(locParticles[loc_j]); |
| 653 | locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locChargedTrackHypothesis, locMatchFOM); |
| 654 | if((locMCThrown == NULL__null) || (locMatchFOM < dMinThrownMatchFOM)) |
| 655 | return false; //not matched |
| 656 | if(((Particle_t)locMCThrown->type) != locParticles[loc_j]->PID()) |
| 657 | return false; //bad PID |
| 658 | } |
| 659 | |
| 660 | //check if parent is correct |
| 661 | auto locParentSearchReactionStep = locReactionStep; |
| 662 | int locParentID = locMCThrown->parentid; |
| 663 | int locSearchStepIndex = loc_i; |
| 664 | |
| 665 | do |
| 666 | { |
| 667 | if(locParentID == -1) |
| 668 | return false; //parent particle is not listed: matched to knock-out particle, is wrong. bail |
| 669 | if(locParentID == 0) |
| 670 | { |
| 671 | //parent id of 0 is directly produced |
| 672 | if(locParentSearchReactionStep->Get_InitialPID() != Gamma) |
| 673 | return false; //was directly (photo-) produced, but not so in combo: bail |
| 674 | break; //good: this is the only "good" exit point of the do-loop |
| 675 | } |
| 676 | |
| 677 | const DMCThrown* locMCThrownParent = locMCThrownMyIDMap[locParentID]; |
| 678 | Particle_t locPID = locMCThrownParent->PID(); |
| 679 | if((locPID == Unknown) || IsResonance(locPID)) |
| 680 | { |
| 681 | //intermediate (unknown or resonance) particle: go to its parent |
| 682 | locParentID = locMCThrownParent->parentid; |
| 683 | continue; |
| 684 | } |
| 685 | if(locPID != locParentSearchReactionStep->Get_InitialPID()) |
| 686 | return false; //the true particle was produced from a different parent: bail |
| 687 | |
| 688 | //OK, we've determined that the particle in question decayed from the correct particle. |
| 689 | //HOWEVER, we need to determine whether the PARENT decayed from the correct particle (and on(back)wards until the production step) |
| 690 | locParentID = locMCThrownParent->parentid; |
| 691 | locSearchStepIndex = DAnalysis::Get_InitialParticleDecayFromIndices(Get_Reaction(), locSearchStepIndex).first; |
| 692 | if(locSearchStepIndex < 0) |
| 693 | { |
| 694 | //DReaction is not the full reaction (i.e. doesn't contain beam (e.g. only pi0 decay)) |
| 695 | if(dExclusiveMatchFlag) |
| 696 | return false; |
| 697 | break; //good |
| 698 | } |
| 699 | locParentSearchReactionStep = Get_Reaction()->Get_ReactionStep(locSearchStepIndex); |
| 700 | } |
| 701 | while(true); |
| 702 | } |
| 703 | } |
| 704 | |
| 705 | return true; //we made it! |
| 706 | } |
| 707 | |
| 708 | DCutAction_TrueCombo::~DCutAction_TrueCombo(void) |
| 709 | { |
| 710 | if(dCutAction_TrueBeamParticle != NULL__null) |
| 711 | delete dCutAction_TrueBeamParticle; |
| 712 | if(dCutAction_ThrownTopology != NULL__null) |
| 713 | delete dCutAction_ThrownTopology; |
| 714 | } |
| 715 | |
| 716 | bool DCutAction_TrueBeamParticle::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 717 | { |
| 718 | vector<const DBeamPhoton*> locBeamPhotons; |
| 719 | locEventLoop->Get(locBeamPhotons, "TAGGEDMCGEN"); |
| 720 | if(locBeamPhotons.empty()) |
| 721 | return false; //true not tagged |
| 722 | |
| 723 | const DKinematicData* locKinematicData = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle_Measured(); |
| 724 | if(locKinematicData == NULL__null) |
| 725 | return false; //initial step is not production step |
| 726 | |
| 727 | const DBeamPhoton* locBeamPhoton = dynamic_cast<const DBeamPhoton*>(locKinematicData); |
| 728 | if(locBeamPhoton == NULL__null) |
| 729 | return false; //dunno how could be possible ... |
| 730 | |
| 731 | double locDeltaT = fabs(locBeamPhoton->time() - locBeamPhotons[0]->time()); |
| 732 | return ((locBeamPhoton->dSystem == locBeamPhotons[0]->dSystem) && (locBeamPhoton->dCounter == locBeamPhotons[0]->dCounter) && (locDeltaT < 1.0)); |
| 733 | } |
| 734 | |
| 735 | bool DCutAction_TruePID::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 736 | { |
| 737 | vector<const DMCThrownMatching*> locMCThrownMatchingVector; |
| 738 | locEventLoop->Get(locMCThrownMatchingVector); |
| 739 | const DMCThrownMatching* locMCThrownMatching = locMCThrownMatchingVector[0]; |
| 740 | |
| 741 | auto locSteps = locParticleCombo->Get_ParticleComboSteps(); |
| 742 | for(size_t loc_i = 0; loc_i < locSteps.size(); ++loc_i) |
| 743 | { |
| 744 | if((dInitialPID != Unknown) && (Get_Reaction()->Get_ReactionStep(loc_i)->Get_InitialPID() != dInitialPID)) |
| 745 | continue; |
| 746 | auto locParticles = locSteps[loc_i]->Get_FinalParticles_Measured(Get_Reaction()->Get_ReactionStep(loc_i), d_AllCharges); |
| 747 | for(size_t loc_j = 0; loc_j < locParticles.size(); ++loc_j) |
| 748 | { |
| 749 | if((locParticles[loc_j]->PID() != dTruePID) && (dTruePID != Unknown)) |
| 750 | continue; |
| 751 | |
| 752 | for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) |
| 753 | { |
| 754 | if(ParticleCharge(dTruePID) == 0) |
| 755 | { |
| 756 | double locMatchFOM = 0.0; |
| 757 | const DNeutralParticleHypothesis* locNeutralParticleHypothesis = static_cast<const DNeutralParticleHypothesis*>(locParticles[loc_i]); |
| 758 | auto locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locNeutralParticleHypothesis, locMatchFOM); |
| 759 | if((locMCThrown == NULL__null) || (locMatchFOM < dMinThrownMatchFOM)) |
| 760 | return false; |
| 761 | if(((Particle_t)locMCThrown->type) != dTruePID) |
| 762 | return false; |
| 763 | } |
| 764 | else |
| 765 | { |
| 766 | double locMatchFOM = 0.0; |
| 767 | const DChargedTrackHypothesis* locChargedTrackHypothesis = static_cast<const DChargedTrackHypothesis*>(locParticles[loc_i]); |
| 768 | auto locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locChargedTrackHypothesis, locMatchFOM); |
| 769 | if((locMCThrown == NULL__null) || (locMatchFOM < dMinThrownMatchFOM)) |
| 770 | return false; |
| 771 | if(((Particle_t)locMCThrown->type) != dTruePID) |
| 772 | return false; |
| 773 | } |
| 774 | } |
| 775 | } |
| 776 | } |
| 777 | return true; |
| 778 | } |
| 779 | |
| 780 | bool DCutAction_AllTruePID::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 781 | { |
| 782 | vector<const DMCThrownMatching*> locMCThrownMatchingVector; |
| 783 | locEventLoop->Get(locMCThrownMatchingVector); |
| 784 | const DMCThrownMatching* locMCThrownMatching = locMCThrownMatchingVector[0]; |
| 785 | const DMCThrown* locMCThrown; |
| 786 | |
| 787 | auto locParticles = locParticleCombo->Get_FinalParticles_Measured(Get_Reaction(), d_AllCharges); |
| 788 | for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) |
| 789 | { |
| 790 | if(ParticleCharge(locParticles[loc_i]->PID()) == 0) |
| 791 | { |
| 792 | double locMatchFOM = 0.0; |
| 793 | const DNeutralParticleHypothesis* locNeutralParticleHypothesis = static_cast<const DNeutralParticleHypothesis*>(locParticles[loc_i]); |
| 794 | locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locNeutralParticleHypothesis, locMatchFOM); |
| 795 | if((locMCThrown == NULL__null) || (locMatchFOM < dMinThrownMatchFOM)) |
| 796 | return false; |
| 797 | if(((Particle_t)locMCThrown->type) != locParticles[loc_i]->PID()) |
| 798 | return false; |
| 799 | } |
| 800 | else |
| 801 | { |
| 802 | double locMatchFOM = 0.0; |
| 803 | const DChargedTrackHypothesis* locChargedTrackHypothesis = static_cast<const DChargedTrackHypothesis*>(locParticles[loc_i]); |
| 804 | locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locChargedTrackHypothesis, locMatchFOM); |
| 805 | if((locMCThrown == NULL__null) || (locMatchFOM < dMinThrownMatchFOM)) |
| 806 | return false; |
| 807 | if(((Particle_t)locMCThrown->type) != locParticles[loc_i]->PID()) |
| 808 | return false; |
| 809 | } |
| 810 | } |
| 811 | return true; |
| 812 | } |
| 813 | |
| 814 | string DCutAction_GoodEventRFBunch::Get_ActionName(void) const |
| 815 | { |
| 816 | ostringstream locStream; |
| 817 | locStream << DAnalysisAction::Get_ActionName() << "_" << dCutIfBadRFBunchFlag; |
| 818 | return locStream.str(); |
| 819 | } |
| 820 | |
| 821 | bool DCutAction_GoodEventRFBunch::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 822 | { |
| 823 | const DEventRFBunch* locEventRFBunch = locParticleCombo->Get_EventRFBunch(); |
| 824 | return (locEventRFBunch->dTime == locEventRFBunch->dTime); |
| 825 | } |
| 826 | |
| 827 | string DCutAction_TransverseMomentum::Get_ActionName(void) const |
| 828 | { |
| 829 | ostringstream locStream; |
| 830 | locStream << DAnalysisAction::Get_ActionName() << "_" << dMaxTransverseMomentum; |
| 831 | return locStream.str(); |
| 832 | } |
| 833 | |
| 834 | bool DCutAction_TransverseMomentum::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 835 | { |
| 836 | auto locParticles = locParticleCombo->Get_FinalParticles_Measured(Get_Reaction(), d_AllCharges); |
| 837 | |
| 838 | DVector3 locTotalMomentum; |
| 839 | for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) |
| 840 | locTotalMomentum += locParticles[loc_i]->momentum(); |
| 841 | |
| 842 | return (dMaxTransverseMomentum >= locTotalMomentum.Perp()); |
| 843 | } |
| 844 | |
| 845 | string DCutAction_TrackHitPattern::Get_ActionName(void) const |
| 846 | { |
| 847 | ostringstream locStream; |
| 848 | locStream << DAnalysisAction::Get_ActionName() << "_" << dMinHitRingsPerCDCSuperlayer << "_" << dMinHitPlanesPerFDCPackage; |
| 849 | return locStream.str(); |
| 850 | } |
| 851 | |
| 852 | bool DCutAction_TrackHitPattern::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 853 | { |
| 854 | auto locParticles = locParticleCombo->Get_FinalParticles_Measured(Get_Reaction(), d_AllCharges); |
| 855 | |
| 856 | const DParticleID* locParticleID = NULL__null; |
| 857 | locEventLoop->GetSingle(locParticleID); |
| 858 | |
| 859 | for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) |
| 860 | { |
| 861 | const DChargedTrackHypothesis* locChargedTrackHypothesis = static_cast<const DChargedTrackHypothesis*>(locParticles[loc_i]); |
| 862 | auto locTrackTimeBased = locChargedTrackHypothesis->Get_TrackTimeBased(); |
| 863 | if(!Cut_TrackHitPattern(locParticleID, locTrackTimeBased)) |
| 864 | return false; |
| 865 | } |
| 866 | |
| 867 | return true; |
| 868 | } |
| 869 | |
| 870 | bool DCutAction_TrackHitPattern::Cut_TrackHitPattern(const DParticleID* locParticleID, const DKinematicData* locTrack) const |
| 871 | { |
| 872 | const DTrackTimeBased* locTrackTimeBased = dynamic_cast<const DTrackTimeBased*>(locTrack); |
| 873 | const DTrackWireBased* locTrackWireBased = dynamic_cast<const DTrackWireBased*>(locTrack); |
| 874 | const DTrackCandidate* locTrackCandidate = dynamic_cast<const DTrackCandidate*>(locTrack); |
| 875 | |
| 876 | map<int, int> locNumHitRingsPerSuperlayer, locNumHitPlanesPerPackage; |
| 877 | if(locTrackTimeBased != NULL__null) |
| 878 | { |
| 879 | locParticleID->Get_CDCNumHitRingsPerSuperlayer(locTrackTimeBased->dCDCRings, locNumHitRingsPerSuperlayer); |
| 880 | locParticleID->Get_FDCNumHitPlanesPerPackage(locTrackTimeBased->dFDCPlanes, locNumHitPlanesPerPackage); |
| 881 | } |
| 882 | else if(locTrackWireBased != NULL__null) |
| 883 | { |
| 884 | locParticleID->Get_CDCNumHitRingsPerSuperlayer(locTrackWireBased->dCDCRings, locNumHitRingsPerSuperlayer); |
| 885 | locParticleID->Get_FDCNumHitPlanesPerPackage(locTrackWireBased->dFDCPlanes, locNumHitPlanesPerPackage); |
| 886 | } |
| 887 | else if(locTrackCandidate != NULL__null) |
| 888 | { |
| 889 | locParticleID->Get_CDCNumHitRingsPerSuperlayer(locTrackCandidate->dCDCRings, locNumHitRingsPerSuperlayer); |
| 890 | locParticleID->Get_FDCNumHitPlanesPerPackage(locTrackCandidate->dFDCPlanes, locNumHitPlanesPerPackage); |
| 891 | } |
| 892 | else |
| 893 | return false; |
| 894 | |
| 895 | //CDC: find inner-most & outer-most superlayers |
| 896 | int locInnermostCDCSuperlayer = 10, locOutermostCDCSuperlayer = 0; |
| 897 | for(auto& locSuperlayerPair : locNumHitRingsPerSuperlayer) |
| 898 | { |
| 899 | if(locSuperlayerPair.second == 0) |
| 900 | continue; //0 hits |
| 901 | if(locSuperlayerPair.first < locInnermostCDCSuperlayer) |
| 902 | locInnermostCDCSuperlayer = locSuperlayerPair.first; |
| 903 | if(locSuperlayerPair.first > locOutermostCDCSuperlayer) |
| 904 | locOutermostCDCSuperlayer = locSuperlayerPair.first; |
| 905 | } |
| 906 | |
| 907 | //CDC: loop again, cutting |
| 908 | for(auto& locSuperlayerPair : locNumHitRingsPerSuperlayer) |
| 909 | { |
| 910 | if(locSuperlayerPair.first == locOutermostCDCSuperlayer) |
| 911 | break; //don't check the last one: track could be leaving |
| 912 | if(locSuperlayerPair.first < locInnermostCDCSuperlayer) |
| 913 | continue; //don't check before the first one: could be detached vertex |
| 914 | if(locSuperlayerPair.second < int(dMinHitRingsPerCDCSuperlayer)) |
| 915 | return false; |
| 916 | } |
| 917 | |
| 918 | //FDC: find inner-most & outer-most superlayers |
| 919 | int locOutermostFDCPlane = 0; |
| 920 | for(auto& locPackagePair : locNumHitPlanesPerPackage) |
| 921 | { |
| 922 | if(locPackagePair.second == 0) |
| 923 | continue; //0 hits |
| 924 | if(locPackagePair.first > locOutermostFDCPlane) |
| 925 | locOutermostFDCPlane = locPackagePair.first; |
| 926 | } |
| 927 | |
| 928 | //FDC: loop again, cutting |
| 929 | for(auto& locPackagePair : locNumHitPlanesPerPackage) |
| 930 | { |
| 931 | if(locPackagePair.first == locOutermostFDCPlane) |
| 932 | break; //don't check the last one: track could be leaving |
| 933 | if(locPackagePair.second == 0) |
| 934 | continue; //0 hits: is ok: could be curling through beamline |
| 935 | if(locPackagePair.second < int(dMinHitPlanesPerFDCPackage)) |
| 936 | return false; |
| 937 | } |
| 938 | |
| 939 | if((locOutermostCDCSuperlayer <= 2) && locNumHitPlanesPerPackage.empty()) |
| 940 | return false; //would have at least expected it to hit the first FDC package //is likely spurious |
| 941 | |
| 942 | return true; |
| 943 | } |
| 944 | |
| 945 | void DCutAction_dEdx::Initialize(JEventLoop* locEventLoop) |
| 946 | { |
| 947 | japp->RootWriteLock(); //ACQUIRE ROOT LOCK!! //I have no idea why this is needed, but without it it crashes. Sigh. |
| 948 | { |
| 949 | if(dCutMap.find(Proton) == dCutMap.end()) |
| 950 | { |
| 951 | dCutMap[Proton].first = new TF1("df_dEdxCut_ProtonLow", "exp(-1.0*[0]*x + [1]) + [2]", 0.0, 12.0); |
| 952 | dCutMap[Proton].first->SetParameters(3.93024, 3.0, 1.0); |
| 953 | dCutMap[Proton].second = new TF1("df_dEdxCut_ProtonHigh", "[0]", 0.0, 12.0); |
| 954 | dCutMap[Proton].second->SetParameter(0, 9999999.9); |
| 955 | } |
| 956 | |
| 957 | if(dCutMap.find(PiPlus) == dCutMap.end()) |
| 958 | { |
| 959 | dCutMap[PiPlus].first = new TF1("df_dEdxCut_PionLow", "[0]", 0.0, 12.0); |
| 960 | dCutMap[PiPlus].first->SetParameter(0, -1.0); |
| 961 | dCutMap[PiPlus].second = new TF1("df_dEdxCut_PionHigh", "exp(-1.0*[0]*x + [1]) + [2]", 0.0, 12.0); |
| 962 | dCutMap[PiPlus].second->SetParameters(6.0, 2.80149, 2.55); |
| 963 | } |
| 964 | } |
| 965 | japp->RootUnLock(); //RELEASE ROOT LOCK!! |
| 966 | } |
| 967 | |
| 968 | bool DCutAction_dEdx::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 969 | { |
| 970 | auto locParticles = locParticleCombo->Get_FinalParticles_Measured(Get_Reaction(), d_Charged); |
| 971 | for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) |
| 972 | { |
| 973 | auto locChargedTrackHypothesis = static_cast<const DChargedTrackHypothesis*>(locParticles[loc_i]); |
| 974 | if(!Cut_dEdx(locChargedTrackHypothesis)) |
| 975 | return false; |
| 976 | } |
| 977 | |
| 978 | return true; |
| 979 | } |
| 980 | |
| 981 | bool DCutAction_dEdx::Cut_dEdx(const DChargedTrackHypothesis* locChargedTrackHypothesis) |
| 982 | { |
| 983 | auto locTrackTimeBased = locChargedTrackHypothesis->Get_TrackTimeBased(); |
| 984 | |
| 985 | Particle_t locPID = locChargedTrackHypothesis->PID(); |
| 986 | if(dCutMap.find(locPID) == dCutMap.end()) |
| 987 | return true; |
| 988 | auto locCutPair = dCutMap[locPID]; |
| 989 | |
| 990 | if(locTrackTimeBased->dNumHitsUsedFordEdx_CDC == 0) |
| 991 | return true; |
| 992 | |
| 993 | auto locP = locTrackTimeBased->momentum().Mag(); |
| 994 | auto locdEdx = locTrackTimeBased->ddEdx_CDC_amp*1.0E6; |
| 995 | |
| 996 | return ((locdEdx >= locCutPair.first->Eval(locP)) && (locdEdx <= locCutPair.second->Eval(locP))); |
| 997 | } |
| 998 | |
| 999 | string DCutAction_BeamEnergy::Get_ActionName(void) const |
| 1000 | { |
| 1001 | ostringstream locStream; |
| 1002 | locStream << DAnalysisAction::Get_ActionName() << "_" << dMinBeamEnergy << "_" << dMaxBeamEnergy; |
| 1003 | return locStream.str(); |
| 1004 | } |
| 1005 | |
| 1006 | bool DCutAction_BeamEnergy::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 1007 | { |
| 1008 | const DKinematicData* locInitParticle = NULL__null; |
| 1009 | if(Get_UseKinFitResultsFlag()) |
| 1010 | locInitParticle = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle(); |
| 1011 | else |
| 1012 | locInitParticle = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle_Measured(); |
| 1013 | if(locInitParticle == NULL__null) |
| 1014 | return false; |
| 1015 | |
| 1016 | double locBeamEnergy = locInitParticle->energy(); |
| 1017 | return ((locBeamEnergy >= dMinBeamEnergy) && (locBeamEnergy <= dMaxBeamEnergy)); |
| 1018 | } |
| 1019 | |
| 1020 | string DCutAction_TrackFCALShowerEOverP::Get_ActionName(void) const |
| 1021 | { |
| 1022 | ostringstream locStream; |
| 1023 | locStream << DAnalysisAction::Get_ActionName() << "_" << dShowerEOverPCut; |
| 1024 | return locStream.str(); |
| 1025 | } |
| 1026 | |
| 1027 | bool DCutAction_TrackFCALShowerEOverP::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 1028 | { |
| 1029 | // For all charged tracks except e+/e-, cuts those with E/p > input value |
| 1030 | // For e+/e-, cuts those with E/p < input value |
| 1031 | // Does not cut tracks without a matching FCAL shower |
| 1032 | |
| 1033 | auto locParticles = Get_UseKinFitResultsFlag() ? locParticleCombo->Get_FinalParticles(Get_Reaction(), false, false, d_Charged) : locParticleCombo->Get_FinalParticles_Measured(Get_Reaction(), d_Charged); |
| 1034 | for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) |
| 1035 | { |
| 1036 | const DChargedTrackHypothesis* locChargedTrackHypothesis = static_cast<const DChargedTrackHypothesis*>(locParticles[loc_i]); |
| 1037 | auto locFCALShowerMatchParams = locChargedTrackHypothesis->Get_FCALShowerMatchParams(); |
| 1038 | if(locFCALShowerMatchParams == NULL__null) |
| 1039 | continue; |
| 1040 | |
| 1041 | Particle_t locPID = locChargedTrackHypothesis->PID(); |
| 1042 | const DFCALShower* locFCALShower = locFCALShowerMatchParams->dFCALShower; |
| 1043 | double locShowerEOverP = locFCALShower->getEnergy()/locChargedTrackHypothesis->momentum().Mag(); |
| 1044 | |
| 1045 | if((locPID == Electron) || (locPID == Positron)) |
| 1046 | { |
| 1047 | if(locShowerEOverP < dShowerEOverPCut) |
| 1048 | return false; |
| 1049 | } |
| 1050 | else if(locShowerEOverP > dShowerEOverPCut) |
| 1051 | return false; |
| 1052 | } |
| 1053 | |
| 1054 | return true; |
| 1055 | } |
| 1056 | |
| 1057 | string DCutAction_TrackShowerEOverP::Get_ActionName(void) const |
| 1058 | { |
| 1059 | ostringstream locStream; |
| 1060 | locStream << DAnalysisAction::Get_ActionName() << "_" << dDetector << "_" << dPID << "_" << dShowerEOverPCut; |
| 1061 | return locStream.str(); |
| 1062 | } |
| 1063 | |
| 1064 | bool DCutAction_TrackShowerEOverP::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 1065 | { |
| 1066 | // For all charged tracks except e+/e-, cuts those with E/p > input value |
| 1067 | // For e+/e-, cuts those with E/p < input value |
| 1068 | // Does not cut tracks without a matching FCAL shower |
| 1069 | |
| 1070 | auto locParticles = Get_UseKinFitResultsFlag() ? locParticleCombo->Get_FinalParticles(Get_Reaction(), false, false, d_Charged) : locParticleCombo->Get_FinalParticles_Measured(Get_Reaction(), d_Charged); |
| 1071 | for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) |
| 1072 | { |
| 1073 | const DChargedTrackHypothesis* locChargedTrackHypothesis = static_cast<const DChargedTrackHypothesis*>(locParticles[loc_i]); |
| 1074 | |
| 1075 | Particle_t locPID = locChargedTrackHypothesis->PID(); |
| 1076 | if(locPID != dPID) |
| 1077 | continue; |
| 1078 | |
| 1079 | double locP = locChargedTrackHypothesis->momentum().Mag(); |
| 1080 | double locShowerEOverP = 0.0; |
| 1081 | if(dDetector == SYS_FCAL) |
| 1082 | { |
| 1083 | auto locFCALShowerMatchParams = locChargedTrackHypothesis->Get_FCALShowerMatchParams(); |
| 1084 | if(locFCALShowerMatchParams == NULL__null) |
| 1085 | continue; |
| 1086 | |
| 1087 | const DFCALShower* locFCALShower = locFCALShowerMatchParams->dFCALShower; |
| 1088 | locShowerEOverP = locFCALShower->getEnergy()/locP; |
| 1089 | } |
| 1090 | else if(dDetector == SYS_BCAL) |
| 1091 | { |
| 1092 | auto locBCALShowerMatchParams = locChargedTrackHypothesis->Get_BCALShowerMatchParams(); |
| 1093 | if(locBCALShowerMatchParams == NULL__null) |
| 1094 | continue; |
| 1095 | |
| 1096 | const DBCALShower* locBCALShower = locBCALShowerMatchParams->dBCALShower; |
| 1097 | locShowerEOverP = locBCALShower->E/locP; |
| 1098 | } |
| 1099 | else |
| 1100 | continue; //what?? |
| 1101 | |
| 1102 | if((locPID == Electron) || (locPID == Positron)) |
| 1103 | { |
| 1104 | if(locShowerEOverP < dShowerEOverPCut) |
| 1105 | return false; |
| 1106 | } |
| 1107 | else if(locShowerEOverP > dShowerEOverPCut) |
| 1108 | return false; |
| 1109 | } |
| 1110 | |
| 1111 | return true; |
| 1112 | } |
| 1113 | |
| 1114 | string DCutAction_PIDDeltaT::Get_ActionName(void) const |
| 1115 | { |
| 1116 | ostringstream locStream; |
| 1117 | locStream << DAnalysisAction::Get_ActionName() << "_" << dPID << "_" << dSystem << "_" << dDeltaTCut; |
| 1118 | return locStream.str(); |
| 1119 | } |
| 1120 | |
| 1121 | bool DCutAction_PIDDeltaT::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 1122 | { |
| 1123 | //if dPID = Unknown, apply cut to all PIDs |
| 1124 | //if dSystem = SYS_NULL, apply cut to all systems |
| 1125 | |
| 1126 | auto locParticles = Get_UseKinFitResultsFlag() ? locParticleCombo->Get_FinalParticles(Get_Reaction(), false, false, d_AllCharges) : locParticleCombo->Get_FinalParticles_Measured(Get_Reaction(), d_AllCharges); |
| 1127 | for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) |
| 1128 | { |
| 1129 | if((dPID != Unknown) && (locParticles[loc_i]->PID() != dPID)) |
| 1130 | continue; |
| 1131 | |
| 1132 | auto locChargedHypo = dynamic_cast<const DChargedTrackHypothesis*>(locParticles[loc_i]); |
| 1133 | if(locChargedHypo != nullptr) |
| 1134 | { |
| 1135 | if((dSystem != SYS_NULL) && (locChargedHypo->t1_detector() != dSystem)) |
| 1136 | continue; |
| 1137 | // double locDeltaT = locChargedHypo->time() - (locChargedHypo->t0() + (locChargedHypo->position().Z() - locVertex.Z())/SPEED_OF_LIGHT); //COMPARE: to old |
| 1138 | double locDeltaT = locChargedHypo->Get_TimeAtPOCAToVertex() - locChargedHypo->t0(); //UNCOMMENT WHEN DONE COMPARING |
| 1139 | if(fabs(locDeltaT) > dDeltaTCut) |
| 1140 | return false; |
| 1141 | continue; |
| 1142 | } |
| 1143 | auto locNeutralHypo = dynamic_cast<const DNeutralParticleHypothesis*>(locParticles[loc_i]); |
| 1144 | if(locNeutralHypo != nullptr) |
| 1145 | { |
| 1146 | if((dSystem != SYS_NULL) && (locNeutralHypo->t1_detector() != dSystem)) |
| 1147 | continue; |
| 1148 | double locDeltaT = locParticles[loc_i]->time() - locNeutralHypo->t0(); |
| 1149 | if(fabs(locDeltaT) > dDeltaTCut) |
| 1150 | return false; |
| 1151 | continue; |
| 1152 | } |
| 1153 | } |
| 1154 | |
| 1155 | return true; |
| 1156 | } |
| 1157 | |
| 1158 | string DCutAction_PIDTimingBeta::Get_ActionName(void) const |
| 1159 | { |
| 1160 | ostringstream locStream; |
| 1161 | locStream << DAnalysisAction::Get_ActionName() << "_" << dPID << "_" << dSystem << "_" << dMinBeta << "_" << dMaxBeta; |
| 1162 | return locStream.str(); |
| 1163 | } |
| 1164 | |
| 1165 | bool DCutAction_PIDTimingBeta::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 1166 | { |
| 1167 | //if dPID = Unknown, apply cut to all PIDs |
| 1168 | //if dSystem = SYS_NULL, apply cut to all systems |
| 1169 | |
| 1170 | auto locParticles = locParticleCombo->Get_FinalParticles_Measured(Get_Reaction(), d_AllCharges); |
| 1171 | for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) |
| 1172 | { |
| 1173 | if((dPID != Unknown) && (locParticles[loc_i]->PID() != dPID)) |
| 1174 | continue; |
| 1175 | |
| 1176 | auto locChargedHypo = dynamic_cast<const DChargedTrackHypothesis*>(locParticles[loc_i]); |
| 1177 | if(locChargedHypo != nullptr) |
| 1178 | { |
| 1179 | if((dSystem != SYS_NULL) && (locChargedHypo->t1_detector() != dSystem)) |
| 1180 | continue; |
| 1181 | double locBeta = locChargedHypo->measuredBeta(); |
| 1182 | if((locBeta < dMinBeta) || (locBeta > dMaxBeta)) |
| 1183 | return false; |
| 1184 | continue; |
| 1185 | } |
| 1186 | auto locNeutralHypo = dynamic_cast<const DNeutralParticleHypothesis*>(locParticles[loc_i]); |
| 1187 | if(locNeutralHypo != nullptr) |
| 1188 | { |
| 1189 | if((dSystem != SYS_NULL) && (locNeutralHypo->t1_detector() != dSystem)) |
| 1190 | continue; |
| 1191 | double locBeta = locNeutralHypo->measuredBeta(); |
| 1192 | if((locBeta < dMinBeta) || (locBeta > dMaxBeta)) |
| 1193 | return false; |
| 1194 | continue; |
| 1195 | } |
| 1196 | } |
| 1197 | |
| 1198 | return true; |
| 1199 | } |
| 1200 | |
| 1201 | string DCutAction_NoPIDHit::Get_ActionName(void) const |
| 1202 | { |
| 1203 | ostringstream locStream; |
| 1204 | locStream << DAnalysisAction::Get_ActionName() << "_" << dPID; |
| 1205 | return locStream.str(); |
| 1206 | } |
| 1207 | |
| 1208 | bool DCutAction_NoPIDHit::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 1209 | { |
| 1210 | //if dPID = Unknown, apply cut to all PIDs |
| 1211 | |
| 1212 | auto locParticles = locParticleCombo->Get_FinalParticles_Measured(Get_Reaction(), d_Charged); |
| 1213 | for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) |
| 1214 | { |
| 1215 | if((dPID != Unknown) && (locParticles[loc_i]->PID() != dPID)) |
| 1216 | continue; |
| 1217 | auto locChargedHypo = static_cast<const DChargedTrackHypothesis*>(locParticles[loc_i]); |
| 1218 | if(locChargedHypo->t1_detector() == SYS_NULL) |
| 1219 | return false; |
| 1220 | } |
| 1221 | |
| 1222 | return true; |
| 1223 | } |
| 1224 | |
| 1225 | string DCutAction_FlightDistance::Get_ActionName(void) const |
| 1226 | { |
| 1227 | ostringstream locStream; |
| 1228 | locStream << DAnalysisAction::Get_ActionName() << "_" << dPID << "_" << dMinFlightDistance; |
| 1229 | return locStream.str(); |
| 1230 | } |
| 1231 | |
| 1232 | bool DCutAction_FlightDistance::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 1233 | { |
| 1234 | //if dPID = Unknown, apply cut to all PIDs |
| 1235 | DKinFitType locKinFitType = Get_Reaction()->Get_KinFitType(); |
| 1236 | |
| 1237 | // for now, require a kinematic fit to make these selections, assuming that the common decay vertex |
| 1238 | // is constrained in the fit |
| 1239 | if(!Get_UseKinFitResultsFlag()) |
| 1240 | return true; |
| 1241 | |
| 1242 | vector<const DReactionVertexInfo*> locVertexInfos; |
| 1243 | locEventLoop->Get(locVertexInfos); |
| 1244 | |
| 1245 | // figure out what the right DReactionVertexInfo is |
| 1246 | const DReactionVertexInfo *locReactionVertexInfo = nullptr; |
| 1247 | for(auto& locVertexInfo : locVertexInfos) { |
| 1248 | auto locVertexReactions = locVertexInfo->Get_Reactions(); |
| 1249 | if(find(locVertexReactions.begin(), locVertexReactions.end(), Get_Reaction()) != locVertexReactions.end()) { |
| 1250 | locReactionVertexInfo = locVertexInfo; |
| 1251 | break; |
| 1252 | } |
| 1253 | } |
| 1254 | |
| 1255 | // check each individual decaying particle (if they exist) |
| 1256 | for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i) |
| 1257 | { |
| 1258 | const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(loc_i); |
| 1259 | //auto locParticles = Get_UseKinFitResultsFlag() ? locParticleComboStep->Get_FinalParticles(Get_Reaction()->Get_ReactionStep(loc_i), false, false, d_AllCharges) : locParticleComboStep->Get_FinalParticles_Measured(Get_Reaction()->Get_ReactionStep(loc_i), d_AllCharges); |
| 1260 | |
| 1261 | if((dPID != Unknown) && (locParticleComboStep->Get_InitialParticle()->PID() != dPID)) |
| 1262 | continue; |
| 1263 | |
| 1264 | // fill info on decaying particles, which is on the particle combo step level |
| 1265 | if(((locParticleComboStep->Get_InitialParticle()!=nullptr) && !Is_FinalStateParticle(locParticleComboStep->Get_InitialParticle()->PID()))) { // decaying particles |
| 1266 | // FOR NOW: we only calculate these displaced vertex quantities if a kinematic fit |
| 1267 | // was performed where the vertex is constrained |
| 1268 | // TODO: Cleverly handle the other cases |
| 1269 | |
| 1270 | // pull out information related to the kinematic fit |
| 1271 | if( !Get_UseKinFitResultsFlag() ) continue; |
| 1272 | if( locReactionVertexInfo == nullptr ) continue; |
| 1273 | |
| 1274 | auto locKinFitParticle = locParticleComboStep->Get_InitialKinFitParticle(); |
| 1275 | |
| 1276 | // extract the info about the vertex, to make sure that the information that we |
| 1277 | // need to calculate displaced vertices is there |
| 1278 | auto locStepVertexInfo = locReactionVertexInfo->Get_StepVertexInfo(loc_i); |
| 1279 | |
| 1280 | auto locPathLength = locKinFitParticle->Get_PathLength(); |
| 1281 | |
| 1282 | if(locPathLength < dMinFlightDistance) |
| 1283 | return false; |
| 1284 | } |
| 1285 | |
| 1286 | } |
| 1287 | |
| 1288 | return true; |
| 1289 | } |
| 1290 | |
| 1291 | |
| 1292 | string DCutAction_FlightSignificance::Get_ActionName(void) const |
| 1293 | { |
| 1294 | ostringstream locStream; |
| 1295 | locStream << DAnalysisAction::Get_ActionName() << "_" << dPID << "_" << dMinFlightSignificance; |
| 1296 | return locStream.str(); |
| 1297 | } |
| 1298 | |
| 1299 | bool DCutAction_FlightSignificance::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 1300 | { |
| 1301 | //if dPID = Unknown, apply cut to all PIDs |
| 1302 | DKinFitType locKinFitType = Get_Reaction()->Get_KinFitType(); |
Value stored to 'locKinFitType' during its initialization is never read | |
| 1303 | |
| 1304 | // for now, require a kinematic fit to make these selections, assuming that the common decay vertex |
| 1305 | // is constrained in the fit |
| 1306 | if(!Get_UseKinFitResultsFlag()) |
| 1307 | return true; |
| 1308 | |
| 1309 | vector<const DReactionVertexInfo*> locVertexInfos; |
| 1310 | locEventLoop->Get(locVertexInfos); |
| 1311 | |
| 1312 | // figure out what the right DReactionVertexInfo is |
| 1313 | const DReactionVertexInfo *locReactionVertexInfo = nullptr; |
| 1314 | for(auto& locVertexInfo : locVertexInfos) { |
| 1315 | auto locVertexReactions = locVertexInfo->Get_Reactions(); |
| 1316 | if(find(locVertexReactions.begin(), locVertexReactions.end(), Get_Reaction()) != locVertexReactions.end()) { |
| 1317 | locReactionVertexInfo = locVertexInfo; |
| 1318 | break; |
| 1319 | } |
| 1320 | } |
| 1321 | |
| 1322 | // check each individual decaying particle (if they exist) |
| 1323 | for(size_t loc_i = 0; loc_i < locParticleCombo->Get_NumParticleComboSteps(); ++loc_i) |
| 1324 | { |
| 1325 | const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(loc_i); |
| 1326 | //auto locParticles = Get_UseKinFitResultsFlag() ? locParticleComboStep->Get_FinalParticles(Get_Reaction()->Get_ReactionStep(loc_i), false, false, d_AllCharges) : locParticleComboStep->Get_FinalParticles_Measured(Get_Reaction()->Get_ReactionStep(loc_i), d_AllCharges); |
| 1327 | |
| 1328 | if((dPID != Unknown) && (locParticleComboStep->Get_InitialParticle()->PID() != dPID)) |
| 1329 | continue; |
| 1330 | |
| 1331 | // fill info on decaying particles, which is on the particle combo step level |
| 1332 | if(((locParticleComboStep->Get_InitialParticle()!=nullptr) && !Is_FinalStateParticle(locParticleComboStep->Get_InitialParticle()->PID()))) { // decaying particles |
| 1333 | // FOR NOW: we only calculate these displaced vertex quantities if a kinematic fit |
| 1334 | // was performed where the vertex is constrained |
| 1335 | // TODO: Cleverly handle the other cases |
| 1336 | |
| 1337 | // pull out information related to the kinematic fit |
| 1338 | if( !Get_UseKinFitResultsFlag() ) continue; |
| 1339 | if( locReactionVertexInfo == nullptr ) continue; |
| 1340 | |
| 1341 | auto locKinFitParticle = locParticleComboStep->Get_InitialKinFitParticle(); |
| 1342 | |
| 1343 | // extract the info about the vertex, to make sure that the information that we |
| 1344 | // need to calculate displaced vertices is there |
| 1345 | auto locStepVertexInfo = locReactionVertexInfo->Get_StepVertexInfo(loc_i); |
| 1346 | |
| 1347 | auto locPathLength = locKinFitParticle->Get_PathLength(); |
| 1348 | auto locPathLengthSigma = locKinFitParticle->Get_PathLengthUncertainty(); |
| 1349 | double locPathLengthSignificance = locPathLength / locPathLengthSigma; |
| 1350 | |
| 1351 | if(locPathLengthSignificance < dMinFlightSignificance) |
| 1352 | return false; |
| 1353 | |
| 1354 | } |
| 1355 | |
| 1356 | } |
| 1357 | |
| 1358 | return true; |
| 1359 | } |
| 1360 | |
| 1361 | |
| 1362 | |
| 1363 | void DCutAction_OneVertexKinFit::Initialize(JEventLoop* locEventLoop) |
| 1364 | { |
| 1365 | // Optional: Useful utility functions. |
| 1366 | locEventLoop->GetSingle(dAnalysisUtilities); |
| 1367 | |
| 1368 | dKinFitUtils = new DKinFitUtils_GlueX(locEventLoop); |
| 1369 | dKinFitter = new DKinFitter(dKinFitUtils); |
| 1370 | dKinFitUtils->Set_UpdateCovarianceMatricesFlag(false); |
| 1371 | |
| 1372 | //CREATE THE HISTOGRAMS |
| 1373 | //Since we are creating histograms, the contents of gDirectory will be modified: must use JANA-wide ROOT lock |
| 1374 | japp->RootWriteLock(); //ACQUIRE ROOT LOCK!! |
| 1375 | { |
| 1376 | //Required: Create a folder in the ROOT output file that will contain all of the output ROOT objects (if any) for this action. |
| 1377 | //If another thread has already created the folder, it just changes to it. |
| 1378 | CreateAndChangeTo_ActionDirectory(); |
| 1379 | |
| 1380 | dHist_ConfidenceLevel = GetOrCreate_Histogram<TH1I>("ConfidenceLevel", "Vertex Kinematic Fit;Confidence Level", 500, 0.0, 1.0); |
| 1381 | dHist_VertexZ = GetOrCreate_Histogram<TH1I>("VertexZ", "Vertex Kinematic Fit;Vertex-Z (cm)", 500, 0.0, 200.0); |
| 1382 | dHist_VertexYVsX = GetOrCreate_Histogram<TH2I>("VertexYVsX", "Vertex Kinematic Fit;Vertex-X (cm);Vertex-Y (cm)", 300, -10.0, 10.0, 300, -10.0, 10.0); |
| 1383 | } |
| 1384 | japp->RootUnLock(); //RELEASE ROOT LOCK!! |
| 1385 | } |
| 1386 | |
| 1387 | void DCutAction_OneVertexKinFit::Run_Update(JEventLoop* locEventLoop) |
| 1388 | { |
| 1389 | locEventLoop->GetSingle(dAnalysisUtilities); |
| 1390 | dKinFitUtils->Set_RunDependent_Data(locEventLoop); |
| 1391 | dKinFitter->Set_RunDependent_Data(locEventLoop); |
| 1392 | } |
| 1393 | |
| 1394 | bool DCutAction_OneVertexKinFit::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) |
| 1395 | { |
| 1396 | //need to call prior to use in each event (cleans up memory allocated from last event) |
| 1397 | //this call invalidates memory from previous fits (but that's OK, we aren't saving them anywhere) |
| 1398 | dKinFitter->Reset_NewEvent(); |
| 1399 | dKinFitUtils->Reset_NewEvent(); |
| 1400 | |
| 1401 | //Get particles for fit (all detected q+) |
| 1402 | auto locParticles = locParticleCombo->Get_FinalParticles_Measured(Get_Reaction(), d_Charged); |
| 1403 | |
| 1404 | //Make DKinFitParticle objects for each one |
| 1405 | deque<shared_ptr<DKinFitParticle>> locKinFitParticles; |
| 1406 | set<shared_ptr<DKinFitParticle>> locKinFitParticleSet; |
| 1407 | for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) |
| 1408 | { |
| 1409 | auto locChargedTrackHypothesis = static_cast<const DChargedTrackHypothesis*>(locParticles[loc_i]); |
| 1410 | auto locKinFitParticle = dKinFitUtils->Make_DetectedParticle(locChargedTrackHypothesis); |
| 1411 | locKinFitParticles.push_back(locKinFitParticle); |
| 1412 | locKinFitParticleSet.insert(locKinFitParticle); |
| 1413 | } |
| 1414 | |
| 1415 | // vertex guess |
| 1416 | TVector3 locVertexGuess = dAnalysisUtilities->Calc_CrudeVertex(locParticles); |
| 1417 | |
| 1418 | // make & set vertex constraint |
| 1419 | auto locVertexConstraint = dKinFitUtils->Make_VertexConstraint(locKinFitParticleSet, {}, locVertexGuess); |
| 1420 | dKinFitter->Add_Constraint(locVertexConstraint); |
| 1421 | |
| 1422 | // PERFORM THE KINEMATIC FIT |
| 1423 | if(!dKinFitter->Fit_Reaction()) |
| 1424 | { |
| 1425 | dKinFitter->Recycle_LastFitMemory(); //RESET MEMORY FROM LAST KINFIT!! //results no longer needed |
| 1426 | return (dMinKinFitCL < 0.0); //fit failed to converge, return false if converge required |
| 1427 | } |
| 1428 | |
| 1429 | // GET THE FIT RESULTS |
| 1430 | double locConfidenceLevel = dKinFitter->Get_ConfidenceLevel(); |
| 1431 | auto locResultVertexConstraint = std::dynamic_pointer_cast<DKinFitConstraint_Vertex>(*dKinFitter->Get_KinFitConstraints().begin()); |
| 1432 | TVector3 locFitVertex = locResultVertexConstraint->Get_CommonVertex(); |
| 1433 | |
| 1434 | //RESET MEMORY FROM LAST KINFIT!! |
| 1435 | dKinFitter->Recycle_LastFitMemory(); //results no longer needed |
| 1436 | |
| 1437 | //FILL HISTOGRAMS |
| 1438 | //Since we are filling histograms local to this action, it will not interfere with other ROOT operations: can use action-wide ROOT lock |
| 1439 | //Note, the mutex is unique to this DReaction + action_string combo: actions of same class with different hists will have a different mutex |
| 1440 | Lock_Action(); //ACQUIRE ROOT LOCK!! |
| 1441 | { |
| 1442 | dHist_ConfidenceLevel->Fill(locConfidenceLevel); |
| 1443 | dHist_VertexZ->Fill(locFitVertex.Z()); |
| 1444 | dHist_VertexYVsX->Fill(locFitVertex.X(), locFitVertex.Y()); |
| 1445 | } |
| 1446 | Unlock_Action(); //RELEASE ROOT LOCK!! |
| 1447 | |
| 1448 | if(locConfidenceLevel < dMinKinFitCL) |
| 1449 | return false; |
| 1450 | if(dMinVertexZ < dMaxVertexZ) //don't cut otherwise |
| 1451 | { |
| 1452 | if((locFitVertex.Z() < dMinVertexZ) || (locFitVertex.Z() > dMaxVertexZ)) |
| 1453 | return false; |
| 1454 | } |
| 1455 | return true; |
| 1456 | } |
| 1457 | |
| 1458 | DCutAction_OneVertexKinFit::~DCutAction_OneVertexKinFit(void) |
| 1459 | { |
| 1460 | if(dKinFitter != NULL__null) |
| 1461 | delete dKinFitter; |
| 1462 | if(dKinFitUtils != NULL__null) |
| 1463 | delete dKinFitUtils; |
| 1464 | } |
| 1465 |