| File: | /home/sdobbs/work/clang/halld_recon/src/libraries/TRACKING/DTrackTimeBased_factory.cc |
| Warning: | line 599, column 7 Called C++ object pointer is null |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
| 1 | // $Id$ | |||
| 2 | // | |||
| 3 | // File: DTrackTimeBased_factory.cc | |||
| 4 | // Created: Thu Sep 4 14:02:44 EDT 2008 | |||
| 5 | // Creator: davidl (on Darwin harriet.jlab.org 8.11.1 i386) | |||
| 6 | // | |||
| 7 | ||||
| 8 | ||||
| 9 | #include <iostream> | |||
| 10 | #include <iomanip> | |||
| 11 | #include <set> | |||
| 12 | #include <mutex> | |||
| 13 | #include <TMath.h> | |||
| 14 | using namespace std; | |||
| 15 | ||||
| 16 | #define TOF_SIGMA0.080 0.080 // TOF resolution in ns | |||
| 17 | ||||
| 18 | #include <TROOT.h> | |||
| 19 | ||||
| 20 | #include "DTrackTimeBased_factory.h" | |||
| 21 | #include <TRACKING/DTrackWireBased.h> | |||
| 22 | #include <TRACKING/DReferenceTrajectory.h> | |||
| 23 | #include <TRACKING/DTrackFitter.h> | |||
| 24 | #include <TRACKING/DTrackHitSelector.h> | |||
| 25 | #include <TRACKING/DMCTrackHit.h> | |||
| 26 | #include <SplitString.h> | |||
| 27 | #include "HDGEOMETRY/DMagneticFieldMapNoField.h" | |||
| 28 | #include <deque> | |||
| 29 | ||||
| 30 | using namespace jana; | |||
| 31 | ||||
| 32 | // Routine for sorting start times | |||
| 33 | bool DTrackTimeBased_T0_cmp(DTrackTimeBased::DStartTime_t a, | |||
| 34 | DTrackTimeBased::DStartTime_t b){ | |||
| 35 | return (a.system>b.system); | |||
| 36 | } | |||
| 37 | ||||
| 38 | bool DTrackTimeBased_cmp(DTrackTimeBased *a,DTrackTimeBased *b){ | |||
| 39 | if (a->candidateid==b->candidateid) return a->mass()<b->mass(); | |||
| 40 | return a->candidateid<b->candidateid; | |||
| 41 | } | |||
| 42 | ||||
| 43 | ||||
| 44 | // count_common_members | |||
| 45 | //------------------ | |||
| 46 | template<typename T> | |||
| 47 | static unsigned int count_common_members(vector<T> &a, vector<T> &b) | |||
| 48 | { | |||
| 49 | unsigned int n=0; | |||
| 50 | for(unsigned int i=0; i<a.size(); i++){ | |||
| 51 | for(unsigned int j=0; j<b.size(); j++){ | |||
| 52 | if(a[i]==b[j])n++; | |||
| 53 | } | |||
| 54 | } | |||
| 55 | ||||
| 56 | return n; | |||
| 57 | } | |||
| 58 | ||||
| 59 | ||||
| 60 | ||||
| 61 | //------------------ | |||
| 62 | // init | |||
| 63 | //------------------ | |||
| 64 | jerror_t DTrackTimeBased_factory::init(void) | |||
| 65 | { | |||
| 66 | fitter = NULL__null; | |||
| 67 | ||||
| 68 | DEBUG_HISTS = false; | |||
| 69 | //DEBUG_HISTS = true; | |||
| 70 | DEBUG_LEVEL = 0; | |||
| 71 | ||||
| 72 | USE_HITS_FROM_WIREBASED_FIT=false; | |||
| 73 | gPARMS->SetDefaultParameter("TRKFIT:USE_HITS_FROM_WIREBASED_FIT", | |||
| 74 | USE_HITS_FROM_WIREBASED_FIT); | |||
| 75 | INSERT_MISSING_HYPOTHESES=true; | |||
| 76 | gPARMS->SetDefaultParameter("TRKFIT:INSERT_MISSING_HYPOTHESES", | |||
| 77 | INSERT_MISSING_HYPOTHESES); | |||
| 78 | ||||
| 79 | gPARMS->SetDefaultParameter("TRKFIT:DEBUG_HISTS",DEBUG_HISTS); | |||
| 80 | gPARMS->SetDefaultParameter("TRKFIT:DEBUG_LEVEL",DEBUG_LEVEL); | |||
| 81 | ||||
| 82 | vector<int> hypotheses; | |||
| 83 | hypotheses.push_back(Positron); | |||
| 84 | hypotheses.push_back(PiPlus); | |||
| 85 | hypotheses.push_back(KPlus); | |||
| 86 | hypotheses.push_back(Proton); | |||
| 87 | hypotheses.push_back(Electron); | |||
| 88 | hypotheses.push_back(PiMinus); | |||
| 89 | hypotheses.push_back(KMinus); | |||
| 90 | hypotheses.push_back(AntiProton); | |||
| 91 | ||||
| 92 | ostringstream locMassStream; | |||
| 93 | for(size_t loc_i = 0; loc_i < hypotheses.size(); ++loc_i) | |||
| 94 | { | |||
| 95 | locMassStream << hypotheses[loc_i]; | |||
| 96 | if(loc_i != (hypotheses.size() - 1)) | |||
| 97 | locMassStream << ","; | |||
| 98 | } | |||
| 99 | ||||
| 100 | string HYPOTHESES = locMassStream.str(); | |||
| 101 | gPARMS->SetDefaultParameter("TRKFIT:HYPOTHESES", HYPOTHESES); | |||
| 102 | ||||
| 103 | // Parse MASS_HYPOTHESES strings to make list of masses to try | |||
| 104 | hypotheses.clear(); | |||
| 105 | SplitString(HYPOTHESES, hypotheses, ","); | |||
| 106 | for(size_t loc_i = 0; loc_i < hypotheses.size(); ++loc_i) | |||
| 107 | { | |||
| 108 | if(ParticleCharge(Particle_t(hypotheses[loc_i])) > 0){ | |||
| 109 | mass_hypotheses_positive.push_back(hypotheses[loc_i]); | |||
| 110 | } | |||
| 111 | else if(ParticleCharge(Particle_t(hypotheses[loc_i])) < 0){ | |||
| 112 | mass_hypotheses_negative.push_back(hypotheses[loc_i]); | |||
| 113 | } | |||
| 114 | } | |||
| 115 | ||||
| 116 | if(mass_hypotheses_positive.empty()){ | |||
| 117 | static once_flag pwarn_flag; | |||
| 118 | call_once(pwarn_flag, [](){ | |||
| 119 | jout << endl; | |||
| 120 | jout << "############# WARNING !! ################ " <<endl; | |||
| 121 | jout << "There are no mass hypotheses for positive tracks!" << endl; | |||
| 122 | jout << "Be SURE this is what you really want!" << endl; | |||
| 123 | jout << "######################################### " <<endl; | |||
| 124 | jout << endl; | |||
| 125 | }); | |||
| 126 | } | |||
| 127 | if(mass_hypotheses_negative.empty()){ | |||
| 128 | static once_flag nwarn_flag; | |||
| 129 | call_once(nwarn_flag, [](){ | |||
| 130 | jout << endl; | |||
| 131 | jout << "############# WARNING !! ################ " <<endl; | |||
| 132 | jout << "There are no mass hypotheses for negative tracks!" << endl; | |||
| 133 | jout << "Be SURE this is what you really want!" << endl; | |||
| 134 | jout << "######################################### " <<endl; | |||
| 135 | jout << endl; | |||
| 136 | }); | |||
| 137 | } | |||
| 138 | ||||
| 139 | mNumHypPlus=mass_hypotheses_positive.size(); | |||
| 140 | mNumHypMinus=mass_hypotheses_negative.size(); | |||
| 141 | ||||
| 142 | // Forces correct particle id (when available) | |||
| 143 | PID_FORCE_TRUTH = false; | |||
| 144 | gPARMS->SetDefaultParameter("TRKFIT:PID_FORCE_TRUTH", PID_FORCE_TRUTH); | |||
| 145 | ||||
| 146 | USE_SC_TIME=true; | |||
| 147 | gPARMS->SetDefaultParameter("TRKFIT:USE_SC_TIME",USE_SC_TIME); | |||
| 148 | ||||
| 149 | USE_TOF_TIME=true; | |||
| 150 | gPARMS->SetDefaultParameter("TRKFIT:USE_TOF_TIME",USE_TOF_TIME); | |||
| 151 | ||||
| 152 | USE_FCAL_TIME=true; | |||
| 153 | gPARMS->SetDefaultParameter("TRKFIT:USE_FCAL_TIME",USE_FCAL_TIME); | |||
| 154 | ||||
| 155 | USE_BCAL_TIME=true; | |||
| 156 | gPARMS->SetDefaultParameter("TRKFIT:USE_BCAL_TIME",USE_BCAL_TIME); | |||
| 157 | ||||
| 158 | return NOERROR; | |||
| 159 | } | |||
| 160 | ||||
| 161 | //------------------ | |||
| 162 | // brun | |||
| 163 | //------------------ | |||
| 164 | jerror_t DTrackTimeBased_factory::brun(jana::JEventLoop *loop, int32_t runnumber) | |||
| 165 | { | |||
| 166 | // Get the geometry | |||
| 167 | DApplication* dapp=dynamic_cast<DApplication*>(loop->GetJApplication()); | |||
| 168 | geom = dapp->GetDGeometry(runnumber); | |||
| 169 | // Check for magnetic field | |||
| 170 | dIsNoFieldFlag = (dynamic_cast<const DMagneticFieldMapNoField*>(dapp->GetBfield(runnumber)) != NULL__null); | |||
| 171 | ||||
| 172 | if(dIsNoFieldFlag){ | |||
| 173 | //Setting this flag makes it so that JANA does not delete the objects in | |||
| 174 | //_data. This factory will manage this memory. | |||
| 175 | //This is all of these pointers are just copied from the "StraightLine" | |||
| 176 | //factory, and should not be re-deleted. | |||
| 177 | SetFactoryFlag(NOT_OBJECT_OWNER); | |||
| 178 | } | |||
| 179 | else{ | |||
| 180 | ClearFactoryFlag(NOT_OBJECT_OWNER); //This factory will create it's own obje | |||
| 181 | } | |||
| 182 | ||||
| 183 | // Get pointer to TrackFitter object that actually fits a track | |||
| 184 | vector<const DTrackFitter *> fitters; | |||
| 185 | loop->Get(fitters); | |||
| 186 | if(fitters.size()<1){ | |||
| 187 | _DBG_std::cerr<<"libraries/TRACKING/DTrackTimeBased_factory.cc" <<":"<<187<<" "<<"Unable to get a DTrackFitter object! NO Charged track fitting will be done!"<<endl; | |||
| 188 | return RESOURCE_UNAVAILABLE; | |||
| 189 | } | |||
| 190 | ||||
| 191 | // Drop the const qualifier from the DTrackFitter pointer (I'm surely going to hell for this!) | |||
| 192 | fitter = const_cast<DTrackFitter*>(fitters[0]); | |||
| 193 | ||||
| 194 | // Warn user if something happened that caused us NOT to get a fitter object pointer | |||
| 195 | if(!fitter){ | |||
| 196 | _DBG_std::cerr<<"libraries/TRACKING/DTrackTimeBased_factory.cc" <<":"<<196<<" "<<"Unable to get a DTrackFitter object! NO Charged track fitting will be done!"<<endl; | |||
| 197 | return RESOURCE_UNAVAILABLE; | |||
| 198 | } | |||
| 199 | ||||
| 200 | // Get the particle ID algorithms | |||
| 201 | vector<const DParticleID *> pid_algorithms; | |||
| 202 | loop->Get(pid_algorithms); | |||
| 203 | if(pid_algorithms.size()<1){ | |||
| 204 | _DBG_std::cerr<<"libraries/TRACKING/DTrackTimeBased_factory.cc" <<":"<<204<<" "<<"Unable to get a DParticleID object! NO PID will be done!"<<endl; | |||
| 205 | return RESOURCE_UNAVAILABLE; | |||
| 206 | } | |||
| 207 | ||||
| 208 | pid_algorithm = pid_algorithms[0]; | |||
| 209 | ||||
| 210 | // Warn user if something happened that caused us NOT to get a pid_algorithm object pointer | |||
| 211 | if(!pid_algorithm){ | |||
| 212 | _DBG_std::cerr<<"libraries/TRACKING/DTrackTimeBased_factory.cc" <<":"<<212<<" "<<"Unable to get a DParticleID object! NO PID will be done!"<<endl; | |||
| 213 | return RESOURCE_UNAVAILABLE; | |||
| 214 | } | |||
| 215 | ||||
| 216 | ||||
| 217 | if(DEBUG_HISTS){ | |||
| 218 | dapp->Lock(); | |||
| 219 | ||||
| 220 | // Histograms may already exist. (Another thread may have created them) | |||
| 221 | // Try and get pointers to the existing ones. | |||
| 222 | fom_chi2_trk = (TH1F*)gROOT(ROOT::GetROOT())->FindObject("fom_chi2_trk"); | |||
| 223 | fom = (TH1F*)gROOT(ROOT::GetROOT())->FindObject("fom"); | |||
| 224 | hitMatchFOM = (TH1F*)gROOT(ROOT::GetROOT())->FindObject("hitMatchFOM"); | |||
| 225 | chi2_trk_mom = (TH2F*)gROOT(ROOT::GetROOT())->FindObject("chi2_trk_mom"); | |||
| 226 | ||||
| 227 | if(!fom_chi2_trk)fom_chi2_trk = new TH1F("fom_chi2_trk","PID FOM: #chi^{2}/Ndf from tracking", 1000, 0.0, 100.0); | |||
| 228 | if(!fom)fom = new TH1F("fom","Combined PID FOM", 1000, 0.0, 1.01); | |||
| 229 | if(!hitMatchFOM)hitMatchFOM = new TH1F("hitMatchFOM","Total Fraction of Hit Matches", 101, 0.0, 1.01); | |||
| 230 | if(!chi2_trk_mom)chi2_trk_mom = new TH2F("chi2_trk_mom","Track #chi^{2}/Ndf versus Kinematic #chi^{2}/Ndf", 1000, 0.0, 100.0, 1000, 0.,100.); | |||
| 231 | ||||
| 232 | ||||
| 233 | Hstart_time=(TH2F*)gROOT(ROOT::GetROOT())->FindObject("Hstart_time"); | |||
| 234 | if (!Hstart_time) Hstart_time=new TH2F("Hstart_time", | |||
| 235 | "vertex time source vs. time", | |||
| 236 | 300,-50,50,9,-0.5,8.5); | |||
| 237 | ||||
| 238 | dapp->Unlock(); | |||
| 239 | ||||
| 240 | } | |||
| 241 | ||||
| 242 | return NOERROR; | |||
| 243 | } | |||
| 244 | ||||
| 245 | //------------------ | |||
| 246 | // evnt | |||
| 247 | //------------------ | |||
| 248 | jerror_t DTrackTimeBased_factory::evnt(JEventLoop *loop, uint64_t eventnumber) | |||
| 249 | { | |||
| 250 | // Save event number to help with debugging | |||
| 251 | myevt=eventnumber; | |||
| 252 | if(!fitter)return NOERROR; | |||
| 253 | ||||
| 254 | if(dIsNoFieldFlag){ | |||
| 255 | //Clear previous objects: | |||
| 256 | //JANA doesn't do it because NOT_OBJECT_OWNER was set | |||
| 257 | //It DID delete them though, in the "StraightLine" factory | |||
| 258 | _data.clear(); | |||
| 259 | ||||
| 260 | vector<const DTrackTimeBased*> locTimeBasedTracks; | |||
| 261 | loop->Get(locTimeBasedTracks, "StraightLine"); | |||
| 262 | for(size_t loc_i = 0; loc_i < locTimeBasedTracks.size(); ++loc_i) | |||
| 263 | _data.push_back(const_cast<DTrackTimeBased*>(locTimeBasedTracks[loc_i])); | |||
| 264 | return NOERROR; | |||
| 265 | } | |||
| 266 | ||||
| 267 | // Get candidates and hits | |||
| 268 | vector<const DTrackWireBased*> tracks; | |||
| 269 | loop->Get(tracks); | |||
| 270 | if (tracks.size()==0) return NOERROR; | |||
| 271 | ||||
| 272 | // get start counter hits | |||
| 273 | vector<const DSCHit*>sc_hits; | |||
| 274 | if (USE_SC_TIME){ | |||
| 275 | loop->Get(sc_hits); | |||
| 276 | } | |||
| 277 | ||||
| 278 | // Get TOF points | |||
| 279 | vector<const DTOFPoint*> tof_points; | |||
| 280 | if (USE_TOF_TIME){ | |||
| 281 | loop->Get(tof_points); | |||
| 282 | } | |||
| 283 | ||||
| 284 | // Get BCAL and FCAL showers | |||
| 285 | vector<const DBCALShower*>bcal_showers; | |||
| 286 | if (USE_BCAL_TIME){ | |||
| 287 | loop->Get(bcal_showers); | |||
| 288 | } | |||
| 289 | vector<const DFCALShower*>fcal_showers; | |||
| 290 | if (USE_FCAL_TIME){ | |||
| 291 | loop->Get(fcal_showers); | |||
| 292 | } | |||
| 293 | ||||
| 294 | vector<const DMCThrown*> mcthrowns; | |||
| 295 | loop->Get(mcthrowns, "FinalState"); | |||
| 296 | ||||
| 297 | // Loop over candidates | |||
| 298 | for(unsigned int i=0; i<tracks.size(); i++){ | |||
| 299 | const DTrackWireBased *track = tracks[i]; | |||
| 300 | ||||
| 301 | unsigned int num=_data.size(); | |||
| 302 | ||||
| 303 | // Create vector of start times from various sources | |||
| 304 | vector<DTrackTimeBased::DStartTime_t>start_times; | |||
| 305 | CreateStartTimeList(track,sc_hits,tof_points,bcal_showers,fcal_showers,start_times); | |||
| 306 | ||||
| 307 | // Fit the track | |||
| 308 | DoFit(track,start_times,loop,track->mass()); | |||
| 309 | ||||
| 310 | //_DBG_<< "eventnumber: " << eventnumber << endl; | |||
| 311 | if (PID_FORCE_TRUTH && _data.size()>num) { | |||
| 312 | // Add figure-of-merit based on difference between thrown and reconstructed momentum | |||
| 313 | // if more than half of the track's hits match MC truth hits and also (charge,mass) | |||
| 314 | // match; add FOM=0 otherwise | |||
| 315 | _data[_data.size()-1]->FOM=GetTruthMatchingFOM(i,_data[_data.size()-1], | |||
| 316 | mcthrowns); | |||
| 317 | } | |||
| 318 | } // loop over track candidates | |||
| 319 | ||||
| 320 | // Filter out duplicate tracks | |||
| 321 | FilterDuplicates(); | |||
| 322 | ||||
| 323 | // Fill in track data for missing hypotheses | |||
| 324 | if (INSERT_MISSING_HYPOTHESES){ | |||
| 325 | InsertMissingHypotheses(loop); | |||
| 326 | } | |||
| 327 | ||||
| 328 | // Set MC Hit-matching information | |||
| 329 | for(size_t loc_i = 0; loc_i < _data.size(); ++loc_i) | |||
| 330 | { | |||
| 331 | if(!mcthrowns.empty()) | |||
| 332 | { | |||
| 333 | double hitFraction; | |||
| 334 | int thrownIndex = GetThrownIndex(mcthrowns, (DKinematicData*)_data[loc_i], hitFraction); | |||
| 335 | _data[loc_i]->dMCThrownMatchMyID = thrownIndex; | |||
| 336 | _data[loc_i]->dNumHitsMatchedToThrown = int(hitFraction * float(Get_NumTrackHits(_data[loc_i])) + 0.01); // + 0.01 so that it rounds down properly | |||
| 337 | } | |||
| 338 | else | |||
| 339 | { | |||
| 340 | _data[loc_i]->dMCThrownMatchMyID = -1; | |||
| 341 | _data[loc_i]->dNumHitsMatchedToThrown = 0; | |||
| 342 | } | |||
| 343 | } | |||
| 344 | ||||
| 345 | return NOERROR; | |||
| 346 | } | |||
| 347 | ||||
| 348 | //------------------ | |||
| 349 | // erun | |||
| 350 | //------------------ | |||
| 351 | jerror_t DTrackTimeBased_factory::erun(void) | |||
| 352 | { | |||
| 353 | return NOERROR; | |||
| 354 | } | |||
| 355 | ||||
| 356 | //------------------ | |||
| 357 | // fini | |||
| 358 | //------------------ | |||
| 359 | jerror_t DTrackTimeBased_factory::fini(void) | |||
| 360 | { | |||
| 361 | return NOERROR; | |||
| 362 | } | |||
| 363 | ||||
| 364 | //------------------ | |||
| 365 | // FilterDuplicates | |||
| 366 | //------------------ | |||
| 367 | void DTrackTimeBased_factory::FilterDuplicates(void) | |||
| 368 | { | |||
| 369 | /// Look through all current DTrackTimeBased objects and remove any | |||
| 370 | /// that have most of their hits in common with another track | |||
| 371 | ||||
| 372 | if(_data.size()==0)return; | |||
| 373 | ||||
| 374 | if(DEBUG_LEVEL>2)_DBG_std::cerr<<"libraries/TRACKING/DTrackTimeBased_factory.cc" <<":"<<374<<" "<<"Looking for clones of time-based tracks ..."<<endl; | |||
| 375 | // We want to remove duplicate tracks corresponding to actual particles, | |||
| 376 | // not just duplicate fitted tracks for certain mass hypotheses -- this | |||
| 377 | // is partly because at a later stage the holes in the list of mass | |||
| 378 | // hypotheses are filled in, thereby spoiling the whole point of this | |||
| 379 | // part of the code! | |||
| 380 | // Keep track of pairs of candididate id's, one for which we want to | |||
| 381 | // keep all the results of fitting with different mass hypotheses, | |||
| 382 | // the other for which we want to delete all the results of fitting. | |||
| 383 | // We need both vectors to take into account potential ambiguities: | |||
| 384 | // for one mass hypothesis starting with one candidate may be "better" | |||
| 385 | // than starting with a second clone candidate, whereas for a second | |||
| 386 | // mass hypothesis, the opposite may be true. | |||
| 387 | vector<unsigned int> candidates_to_keep; | |||
| 388 | vector<unsigned int> candidates_to_delete; | |||
| 389 | for(unsigned int i=0; i<_data.size()-1; i++){ | |||
| 390 | DTrackTimeBased *dtrack1 = _data[i]; | |||
| 391 | ||||
| 392 | vector<const DCDCTrackHit*> cdchits1; | |||
| 393 | vector<const DFDCPseudo*> fdchits1; | |||
| 394 | dtrack1->Get(cdchits1); | |||
| 395 | dtrack1->Get(fdchits1); | |||
| 396 | // Total number of hits in this candidate | |||
| 397 | unsigned int num_cdc1=cdchits1.size(); | |||
| 398 | unsigned int num_fdc1=fdchits1.size(); | |||
| 399 | unsigned int total1 = num_cdc1+num_fdc1; | |||
| 400 | ||||
| 401 | JObject::oid_t cand1=dtrack1->candidateid; | |||
| 402 | for(unsigned int j=i+1; j<_data.size(); j++){ | |||
| 403 | DTrackTimeBased *dtrack2 = _data[j]; | |||
| 404 | if (dtrack2->candidateid==cand1) continue; | |||
| 405 | ||||
| 406 | vector<const DCDCTrackHit*> cdchits2; | |||
| 407 | vector<const DFDCPseudo*> fdchits2; | |||
| 408 | dtrack2->Get(cdchits2); | |||
| 409 | dtrack2->Get(fdchits2); | |||
| 410 | ||||
| 411 | // Total number of hits in this candidate | |||
| 412 | unsigned int num_cdc2=cdchits2.size(); | |||
| 413 | unsigned int num_fdc2=fdchits2.size(); | |||
| 414 | unsigned int total2 = num_cdc2+num_fdc2; | |||
| 415 | ||||
| 416 | // Count number of cdc and fdc hits in common | |||
| 417 | unsigned int Ncdc = count_common_members(cdchits1, cdchits2); | |||
| 418 | unsigned int Nfdc = count_common_members(fdchits1, fdchits2); | |||
| 419 | ||||
| 420 | if(DEBUG_LEVEL>3){ | |||
| 421 | _DBG_std::cerr<<"libraries/TRACKING/DTrackTimeBased_factory.cc" <<":"<<421<<" "<<"cand1:"<<cand1<<" cand2:"<<dtrack2->candidateid<<endl; | |||
| 422 | _DBG_std::cerr<<"libraries/TRACKING/DTrackTimeBased_factory.cc" <<":"<<422<<" "<<" Ncdc="<<Ncdc<<" num_cdc1="<<num_cdc1<<" num_cdc2="<<num_cdc2<<endl; | |||
| 423 | _DBG_std::cerr<<"libraries/TRACKING/DTrackTimeBased_factory.cc" <<":"<<423<<" "<<" Nfdc="<<Nfdc<<" num_fdc1="<<num_fdc1<<" num_fdc2="<<num_fdc2<<endl; | |||
| 424 | } | |||
| 425 | unsigned int total = Ncdc + Nfdc; | |||
| 426 | // If the tracks share at most one hit, consider them | |||
| 427 | // to be separate tracks | |||
| 428 | if (total<=1) continue; | |||
| 429 | ||||
| 430 | // Deal with the case where there are cdc hits in | |||
| 431 | // common between the tracks but there were no fdc | |||
| 432 | // hits used in one of the tracks. | |||
| 433 | if (Ncdc>0 && (num_fdc1>0 || num_fdc2>0) | |||
| 434 | && (num_fdc1*num_fdc2)==0) continue; | |||
| 435 | ||||
| 436 | // Deal with the case where there are fdc hits in | |||
| 437 | // common between the tracks but no cdc hits used in | |||
| 438 | // one of the tracks. | |||
| 439 | if (Nfdc>0 && (num_cdc1>0 || num_cdc2>0) | |||
| 440 | && (num_cdc1*num_cdc2)==0) continue; | |||
| 441 | ||||
| 442 | // Look for tracks with many common hits in the CDC | |||
| 443 | if (num_cdc1>0 && num_cdc2>0){ | |||
| 444 | if (double(Ncdc)/double(num_cdc1)<0.9) continue; | |||
| 445 | if (double(Ncdc)/double(num_cdc2)<0.9) continue; | |||
| 446 | } | |||
| 447 | // Look for tracks with many common hits in the FDC | |||
| 448 | if (num_fdc1>0 && num_fdc2>0){ | |||
| 449 | if (double(Nfdc)/double(num_fdc1)<0.9) continue; | |||
| 450 | if (double(Nfdc)/double(num_fdc2)<0.9) continue; | |||
| 451 | } | |||
| 452 | ||||
| 453 | if(total1<total2){ | |||
| 454 | candidates_to_delete.push_back(cand1); | |||
| 455 | candidates_to_keep.push_back(dtrack2->candidateid); | |||
| 456 | }else if(total2<total1){ | |||
| 457 | candidates_to_delete.push_back(dtrack2->candidateid); | |||
| 458 | candidates_to_keep.push_back(cand1); | |||
| 459 | }else if(dtrack1->FOM > dtrack2->FOM){ | |||
| 460 | candidates_to_delete.push_back(dtrack2->candidateid); | |||
| 461 | candidates_to_keep.push_back(cand1); | |||
| 462 | }else{ | |||
| 463 | candidates_to_delete.push_back(cand1); | |||
| 464 | candidates_to_keep.push_back(dtrack2->candidateid); | |||
| 465 | } | |||
| 466 | } | |||
| 467 | } | |||
| 468 | ||||
| 469 | if(DEBUG_LEVEL>2)_DBG_std::cerr<<"libraries/TRACKING/DTrackTimeBased_factory.cc" <<":"<<469<<" "<<"Found "<<candidates_to_delete.size()<<" time-based clones"<<endl; | |||
| 470 | ||||
| 471 | ||||
| 472 | // Return now if we're keeping everyone | |||
| 473 | if(candidates_to_delete.size()==0)return; | |||
| 474 | ||||
| 475 | // Deal with the ambiguity problem mentioned above | |||
| 476 | for (unsigned int i=0;i<candidates_to_keep.size();i++){ | |||
| 477 | for (unsigned int j=0;j<candidates_to_delete.size();j++){ | |||
| 478 | if (candidates_to_keep[i]==candidates_to_delete[j]){ | |||
| 479 | candidates_to_delete.erase(candidates_to_delete.begin()+j); | |||
| 480 | break; | |||
| 481 | } | |||
| 482 | } | |||
| 483 | ||||
| 484 | } | |||
| 485 | ||||
| 486 | // Copy pointers that we want to keep to a new container and delete | |||
| 487 | // the clone objects | |||
| 488 | vector<DTrackTimeBased*> new_data; | |||
| 489 | sort(_data.begin(),_data.end(),DTrackTimeBased_cmp); | |||
| 490 | for (unsigned int i=0;i<_data.size();i++){ | |||
| 491 | bool keep_track=true; | |||
| 492 | for (unsigned int j=0;j<candidates_to_delete.size();j++){ | |||
| 493 | if (_data[i]->candidateid==candidates_to_delete[j]){ | |||
| 494 | keep_track=false; | |||
| 495 | if(DEBUG_LEVEL>1){ | |||
| 496 | _DBG_std::cerr<<"libraries/TRACKING/DTrackTimeBased_factory.cc" <<":"<<496<<" "<<"Deleting clone time-based fitted result "<<i | |||
| 497 | << " in event " << myevt << endl; | |||
| 498 | } | |||
| 499 | break; | |||
| 500 | } | |||
| 501 | } | |||
| 502 | if (keep_track){ | |||
| 503 | new_data.push_back(_data[i]); | |||
| 504 | } | |||
| 505 | else delete _data[i]; | |||
| 506 | } | |||
| 507 | _data = new_data; | |||
| 508 | } | |||
| 509 | ||||
| 510 | // Returns a FOM based on difference between thrown and reconstructed momentum if track matches MC truth information, | |||
| 511 | // returns a FOM=0 otherwise; | |||
| 512 | // a match requires identical masses and charges, and that more than half of the track's hits match the truth hits | |||
| 513 | double DTrackTimeBased_factory::GetTruthMatchingFOM(int trackIndex,DTrackTimeBased *track,vector<const DMCThrown*>mcthrowns) { | |||
| 514 | bool match=false; | |||
| 515 | ||||
| 516 | DLorentzVector fourMom = track->lorentzMomentum(); | |||
| 517 | //DLorentzVector gen_fourMom[mcthrowns.size()]; | |||
| 518 | vector<DLorentzVector> gen_fourMom(mcthrowns.size()); | |||
| 519 | for(unsigned int i=0; i<mcthrowns.size(); i++){ | |||
| ||||
| 520 | gen_fourMom[i] = mcthrowns[i]->lorentzMomentum(); | |||
| 521 | } | |||
| 522 | ||||
| 523 | // Get info for thrown track | |||
| 524 | double f = 0.; | |||
| 525 | int thrownIndex = GetThrownIndex(mcthrowns, track, f); | |||
| 526 | if(thrownIndex<=0 || f<=0.5) return 0.; | |||
| 527 | ||||
| 528 | double delta_pt_over_pt = (fourMom.Pt()-gen_fourMom[thrownIndex-1].Pt())/gen_fourMom[thrownIndex-1].Pt(); | |||
| 529 | double delta_theta = (fourMom.Theta()-gen_fourMom[thrownIndex-1].Theta())*1000.0; // in milliradians | |||
| 530 | double delta_phi = (fourMom.Phi()-gen_fourMom[thrownIndex-1].Phi())*1000.0; // in milliradians | |||
| 531 | double chisq = pow(delta_pt_over_pt/0.04, 2.0) + pow(delta_theta/20.0, 2.0) + pow(delta_phi/20.0, 2.0); | |||
| 532 | ||||
| 533 | if (fabs(track->mass()-mcthrowns[thrownIndex-1]->mass())<0.01 && track->charge()==mcthrowns[thrownIndex-1]->charge()) | |||
| 534 | match = true; | |||
| 535 | ||||
| 536 | double trk_chi2=track->chisq; | |||
| 537 | unsigned int ndof=track->Ndof; | |||
| 538 | ||||
| 539 | if(DEBUG_HISTS&&match){ | |||
| 540 | fom_chi2_trk->Fill(track->chisq); | |||
| 541 | chi2_trk_mom->Fill(chisq/3.,trk_chi2/ndof); | |||
| 542 | fom->Fill(TMath::Prob(chisq,3)); | |||
| 543 | } | |||
| 544 | ||||
| 545 | /*_DBG_ << "f: " << f << endl; | |||
| 546 | _DBG_ << "trk_chi2: " << trk_chi2 << endl; | |||
| 547 | _DBG_ << "ndof: " << ndof << endl; | |||
| 548 | _DBG_ << "throwncharge: " << mcthrowns[thrownIndex-1]->charge() << endl; | |||
| 549 | _DBG_ << "trackcharge: " << track->charge() << endl; | |||
| 550 | _DBG_ << "chargediff: " << fabs(track->charge()-mcthrowns[thrownIndex-1]->charge()) << endl; | |||
| 551 | _DBG_ << "thrownmass: " << mcthrowns[thrownIndex-1]->mass() << endl; | |||
| 552 | _DBG_ << "trackmass: " << track->mass() << endl; | |||
| 553 | _DBG_ << "massdiff: " << fabs(track->mass()-mcthrowns[thrownIndex-1]->mass()) << endl; | |||
| 554 | _DBG_ << "chisq: " << chisq << endl; | |||
| 555 | _DBG_ << "match?: " << match << endl; | |||
| 556 | _DBG_ << "thrownIndex: " << thrownIndex << " trackIndex: " << trackIndex << endl; | |||
| 557 | _DBG_<< "track " << setprecision(4) << "Px: " << fourMom.Px() << " Py: " << fourMom.Py() << " Pz: " << fourMom.Pz() << " E: " << fourMom.E() << " M: " << fourMom.M() << " pt: " << fourMom.Pt() << " theta: " << fourMom.Theta() << " phi: " << fourMom.Phi() << endl; | |||
| 558 | _DBG_<< "thrown " << setprecision(4) << "Px: " << gen_fourMom[thrownIndex-1].Px() << " Py: " << gen_fourMom[thrownIndex-1].Py() << " Pz: " << gen_fourMom[thrownIndex-1].Pz() << " E: " << gen_fourMom[thrownIndex-1].E() << " M: " << gen_fourMom[thrownIndex-1].M() << " pt: " << gen_fourMom[thrownIndex-1].Pt() << " theta: " << gen_fourMom[thrownIndex-1].Theta() << " phi: " << gen_fourMom[thrownIndex-1].Phi() << endl;*/ | |||
| 559 | ||||
| 560 | return (match) ? TMath::Prob(chisq,3) : 0.0; | |||
| 561 | } | |||
| 562 | ||||
| 563 | //------------------ | |||
| 564 | // GetThrownIndex | |||
| 565 | //------------------ | |||
| 566 | int DTrackTimeBased_factory::GetThrownIndex(vector<const DMCThrown*>& locMCThrowns, const DKinematicData *kd, double &f) | |||
| 567 | { | |||
| 568 | // The DKinematicData object should be a DTrackCandidate, DTrackWireBased, or DParticle which | |||
| 569 | // has associated objects for the hits | |||
| 570 | vector<const DCDCTrackHit*> cdctrackhits; | |||
| 571 | kd->Get(cdctrackhits); | |||
| 572 | vector<const DFDCPseudo*> fdcpseudos; | |||
| 573 | kd->Get(fdcpseudos); | |||
| 574 | ||||
| 575 | int locTotalNumHits = cdctrackhits.size() + fdcpseudos.size(); | |||
| 576 | if(locTotalNumHits == 0) | |||
| 577 | { | |||
| 578 | f = 0; | |||
| 579 | return -1; | |||
| 580 | } | |||
| 581 | ||||
| 582 | // The track number is buried in the truth hit objects of type DMCTrackHit. These should be | |||
| 583 | // associated objects for the individual hit objects. We need to loop through them and | |||
| 584 | // keep track of how many hits for each track number we find | |||
| 585 | ||||
| 586 | map<int, int> locHitMatches; //first int is MC my id, second is num hits | |||
| 587 | for(size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i) | |||
| 588 | { | |||
| 589 | if(fabs(locMCThrowns[loc_i]->charge()) > 0.9) | |||
| 590 | locHitMatches[locMCThrowns[loc_i]->myid] = 0; | |||
| 591 | } | |||
| 592 | ||||
| 593 | // CDC hits | |||
| 594 | for(size_t loc_i = 0; loc_i < cdctrackhits.size(); ++loc_i) | |||
| 595 | { | |||
| 596 | const DCDCHit* locCDCHit = NULL__null; | |||
| 597 | cdctrackhits[loc_i]->GetSingle(locCDCHit); | |||
| 598 | vector<const DCDCHit*> locTruthCDCHits; | |||
| 599 | locCDCHit->Get(locTruthCDCHits); | |||
| ||||
| 600 | if(locTruthCDCHits.empty()) continue; | |||
| 601 | ||||
| 602 | int itrack = locTruthCDCHits[0]->itrack; | |||
| 603 | if(locHitMatches.find(itrack) == locHitMatches.end()) | |||
| 604 | continue; | |||
| 605 | ++locHitMatches[itrack]; | |||
| 606 | } | |||
| 607 | ||||
| 608 | // FDC hits | |||
| 609 | for(size_t loc_i = 0; loc_i < fdcpseudos.size(); ++loc_i) | |||
| 610 | { | |||
| 611 | vector<const DMCTrackHit*> mctrackhits; | |||
| 612 | fdcpseudos[loc_i]->Get(mctrackhits); | |||
| 613 | if(mctrackhits.empty()) | |||
| 614 | continue; | |||
| 615 | if(!mctrackhits[0]->primary) | |||
| 616 | continue; | |||
| 617 | ||||
| 618 | int itrack = mctrackhits[0]->itrack; | |||
| 619 | if(locHitMatches.find(itrack) == locHitMatches.end()) | |||
| 620 | continue; | |||
| 621 | ++locHitMatches[itrack]; | |||
| 622 | } | |||
| 623 | ||||
| 624 | // Find DMCThrown::myid with most wires hit | |||
| 625 | map<int, int>::iterator locIterator = locHitMatches.begin(); | |||
| 626 | int locBestMyID = -1; | |||
| 627 | int locBestNumHits = 0; | |||
| 628 | for(; locIterator != locHitMatches.end(); ++locIterator) | |||
| 629 | { | |||
| 630 | if(locIterator->second <= locBestNumHits) | |||
| 631 | continue; | |||
| 632 | locBestNumHits = locIterator->second; | |||
| 633 | locBestMyID = locIterator->first; | |||
| 634 | } | |||
| 635 | ||||
| 636 | // total fraction of reconstructed hits that match truth hits | |||
| 637 | f = 1.0*locBestNumHits/locTotalNumHits; | |||
| 638 | if(DEBUG_HISTS)hitMatchFOM->Fill(f); | |||
| 639 | ||||
| 640 | return locBestMyID; | |||
| 641 | } | |||
| 642 | ||||
| 643 | ||||
| 644 | // Create a list of start (vertex) times from various sources, ordered by | |||
| 645 | // uncertainty. | |||
| 646 | void DTrackTimeBased_factory | |||
| 647 | ::CreateStartTimeList(const DTrackWireBased *track, | |||
| 648 | vector<const DSCHit*>&sc_hits, | |||
| 649 | vector<const DTOFPoint*>&tof_points, | |||
| 650 | vector<const DBCALShower*>&bcal_showers, | |||
| 651 | vector<const DFCALShower*>&fcal_showers, | |||
| 652 | vector<DTrackTimeBased::DStartTime_t>&start_times){ | |||
| 653 | DTrackTimeBased::DStartTime_t start_time; | |||
| 654 | ||||
| 655 | // Match to the start counter and the outer detectors | |||
| 656 | double locStartTimeVariance = 0.0; | |||
| 657 | double track_t0=track->t0(); | |||
| 658 | double locStartTime = track_t0; // initial guess from tracking | |||
| 659 | ||||
| 660 | // Get start time estimate from Start Counter | |||
| 661 | if (pid_algorithm->Get_StartTime(track->extrapolations.at(SYS_START),sc_hits,locStartTime)){ | |||
| 662 | start_time.t0=locStartTime; | |||
| 663 | // start_time.t0_sigma=sqrt(locTimeVariance); //uncomment when ready | |||
| 664 | start_time.t0_sigma=sqrt(locStartTimeVariance); | |||
| 665 | start_time.system=SYS_START; | |||
| 666 | start_times.push_back(start_time); | |||
| 667 | } | |||
| 668 | // Get start time estimate from TOF | |||
| 669 | locStartTime = track_t0; // initial guess from tracking | |||
| 670 | if (pid_algorithm->Get_StartTime(track->extrapolations.at(SYS_TOF),tof_points,locStartTime)){ | |||
| 671 | // Fill in the start time vector | |||
| 672 | start_time.t0=locStartTime; | |||
| 673 | start_time.t0_sigma=sqrt(locStartTimeVariance); | |||
| 674 | // start_time.t0_sigma=sqrt(locTimeVariance); //uncomment when ready | |||
| 675 | start_time.system=SYS_TOF; | |||
| 676 | start_times.push_back(start_time); | |||
| 677 | } | |||
| 678 | // Get start time estimate from FCAL | |||
| 679 | locStartTime = track_t0; // initial guess from tracking | |||
| 680 | if (pid_algorithm->Get_StartTime(track->extrapolations.at(SYS_FCAL),fcal_showers,locStartTime)){ | |||
| 681 | // Fill in the start time vector | |||
| 682 | start_time.t0=locStartTime; | |||
| 683 | start_time.t0_sigma=sqrt(locStartTimeVariance); | |||
| 684 | // start_time.t0_sigma=sqrt(locTimeVariance); //uncomment when ready | |||
| 685 | start_time.system=SYS_FCAL; | |||
| 686 | start_times.push_back(start_time); | |||
| 687 | } | |||
| 688 | // Get start time estimate from BCAL | |||
| 689 | locStartTime=track_t0; | |||
| 690 | if (pid_algorithm->Get_StartTime(track->extrapolations.at(SYS_BCAL),bcal_showers,locStartTime)){ | |||
| 691 | // Fill in the start time vector | |||
| 692 | start_time.t0=locStartTime; | |||
| 693 | start_time.t0_sigma=0.5; | |||
| 694 | // start_time.t0_sigma=sqrt(locTimeVariance); //uncomment when ready | |||
| 695 | start_time.system=SYS_BCAL; | |||
| 696 | start_times.push_back(start_time); | |||
| 697 | } | |||
| 698 | // Add the t0 estimate from the tracking | |||
| 699 | start_time.t0=track_t0; | |||
| 700 | start_time.t0_sigma=5.; | |||
| 701 | start_time.system=track->t0_detector(); | |||
| 702 | start_times.push_back(start_time); | |||
| 703 | ||||
| 704 | // Set t0 for the fit to the first entry in the list. Usually this will be | |||
| 705 | // from the start counter. | |||
| 706 | mStartTime=start_times[0].t0; | |||
| 707 | mStartDetector=start_times[0].system; | |||
| 708 | ||||
| 709 | // _DBG_ << mStartDetector << " " << mStartTime << endl; | |||
| 710 | ||||
| 711 | } | |||
| 712 | ||||
| 713 | // Create a list of start times and do the fit for a particular mass hypothesis | |||
| 714 | bool DTrackTimeBased_factory::DoFit(const DTrackWireBased *track, | |||
| 715 | vector<DTrackTimeBased::DStartTime_t>&start_times, | |||
| 716 | JEventLoop *loop, | |||
| 717 | double mass){ | |||
| 718 | if(DEBUG_LEVEL>1){_DBG__std::cerr<<"libraries/TRACKING/DTrackTimeBased_factory.cc" <<":"<<718<<std::endl;_DBG_std::cerr<<"libraries/TRACKING/DTrackTimeBased_factory.cc" <<":"<<718<<" "<<"---- Starting time based fit with mass: "<<mass<<endl;} | |||
| 719 | // Get the hits from the wire-based track | |||
| 720 | vector<const DFDCPseudo*>myfdchits; | |||
| 721 | track->GetT(myfdchits); | |||
| 722 | vector<const DCDCTrackHit *>mycdchits; | |||
| 723 | track->GetT(mycdchits); | |||
| 724 | ||||
| 725 | // Do the fit | |||
| 726 | DTrackFitter::fit_status_t status = DTrackFitter::kFitNotDone; | |||
| 727 | if (USE_HITS_FROM_WIREBASED_FIT) { | |||
| 728 | fitter->Reset(); | |||
| 729 | fitter->SetFitType(DTrackFitter::kTimeBased); | |||
| 730 | ||||
| 731 | fitter->AddHits(myfdchits); | |||
| 732 | fitter->AddHits(mycdchits); | |||
| 733 | ||||
| 734 | status=fitter->FitTrack(track->position(),track->momentum(), | |||
| 735 | track->charge(),mass,mStartTime,mStartDetector); | |||
| 736 | } | |||
| 737 | else{ | |||
| 738 | fitter->Reset(); | |||
| 739 | fitter->SetFitType(DTrackFitter::kTimeBased); | |||
| 740 | status = fitter->FindHitsAndFitTrack(*track, track->extrapolations,loop, | |||
| 741 | mass, | |||
| 742 | mycdchits.size()+2*myfdchits.size(), | |||
| 743 | mStartTime,mStartDetector); | |||
| 744 | ||||
| 745 | // If the status is kFitNotDone, then not enough hits were attached to this | |||
| 746 | // track using the hit-gathering algorithm. In this case get the hits | |||
| 747 | // from the wire-based track | |||
| 748 | if (status==DTrackFitter::kFitNotDone){ | |||
| 749 | //_DBG_ << " Using wire-based hits " << endl; | |||
| 750 | fitter->Reset(); | |||
| 751 | fitter->SetFitType(DTrackFitter::kTimeBased); | |||
| 752 | fitter->AddHits(myfdchits); | |||
| 753 | fitter->AddHits(mycdchits); | |||
| 754 | ||||
| 755 | status=fitter->FitTrack(track->position(),track->momentum(), | |||
| 756 | track->charge(),mass,mStartTime,mStartDetector); | |||
| 757 | } | |||
| 758 | ||||
| 759 | } | |||
| 760 | ||||
| 761 | // if the fit returns chisq=-1, something went terribly wrong. We may still | |||
| 762 | // have a usable track from the wire-based pass. In this case set | |||
| 763 | // kFitNoImprovement so we can save the wire-based results. | |||
| 764 | if (fitter->GetChisq()<0){ | |||
| 765 | status=DTrackFitter::kFitNoImprovement; | |||
| 766 | } | |||
| 767 | ||||
| 768 | // In the transition region between the CDC and the FDC where the track | |||
| 769 | // contains both CDC and FDC hits, sometimes too many hits are discarded in | |||
| 770 | // the time-based phase and the time-based fit result does not improve on the | |||
| 771 | // wire-based fit result. In this case set the status word to | |||
| 772 | // kFitNoImprovement and copy the wire-based parameters into the time-based | |||
| 773 | // class. | |||
| 774 | if (myfdchits.size()>3 && mycdchits.size()>3){ | |||
| 775 | unsigned int ndof=fitter->GetNdof(); | |||
| 776 | if (TMath::Prob(track->chisq,track->Ndof)> | |||
| 777 | TMath::Prob(fitter->GetChisq(),ndof)&&ndof<5) | |||
| 778 | status=DTrackFitter::kFitNoImprovement; | |||
| 779 | } | |||
| 780 | ||||
| 781 | // Check the status value from the fit | |||
| 782 | switch(status){ | |||
| 783 | case DTrackFitter::kFitNotDone: | |||
| 784 | //_DBG_<<"Fitter returned kFitNotDone. This should never happen!!"<<endl; | |||
| 785 | case DTrackFitter::kFitFailed: | |||
| 786 | break; | |||
| 787 | case DTrackFitter::kFitNoImprovement: | |||
| 788 | { | |||
| 789 | // Create a new time-based track object | |||
| 790 | DTrackTimeBased *timebased_track = new DTrackTimeBased(); | |||
| 791 | *static_cast<DTrackingData*>(timebased_track) = *static_cast<const DTrackingData*>(track); | |||
| 792 | ||||
| 793 | timebased_track->chisq = track->chisq; | |||
| 794 | timebased_track->Ndof = track->Ndof; | |||
| 795 | timebased_track->pulls = track->pulls; | |||
| 796 | timebased_track->extrapolations = track->extrapolations; | |||
| 797 | timebased_track->trackid = track->id; | |||
| 798 | timebased_track->candidateid=track->candidateid; | |||
| 799 | timebased_track->FOM=track->FOM; | |||
| 800 | timebased_track->flags=DTrackTimeBased::FLAG__USED_WIREBASED_FIT; | |||
| 801 | ||||
| 802 | // add the list of start times | |||
| 803 | timebased_track->start_times.assign(start_times.begin(), | |||
| 804 | start_times.end()); | |||
| 805 | ||||
| 806 | for(unsigned int m=0; m<myfdchits.size(); m++) | |||
| 807 | timebased_track->AddAssociatedObject(myfdchits[m]); | |||
| 808 | for(unsigned int m=0; m<mycdchits.size(); m++) | |||
| 809 | timebased_track->AddAssociatedObject(mycdchits[m]); | |||
| 810 | ||||
| 811 | timebased_track->measured_cdc_hits_on_track = mycdchits.size(); | |||
| 812 | timebased_track->measured_fdc_hits_on_track = myfdchits.size(); | |||
| 813 | ||||
| 814 | // dEdx | |||
| 815 | double locdEdx_FDC, locdx_FDC, locdEdx_CDC, locdEdx_CDC_amp; | |||
| 816 | double locdx_CDC_amp,locdx_CDC; | |||
| 817 | unsigned int locNumHitsUsedFordEdx_FDC, locNumHitsUsedFordEdx_CDC; | |||
| 818 | pid_algorithm->CalcDCdEdx(timebased_track, locdEdx_FDC, locdx_FDC, locdEdx_CDC, locdEdx_CDC_amp, locdx_CDC, locdx_CDC_amp,locNumHitsUsedFordEdx_FDC, locNumHitsUsedFordEdx_CDC); | |||
| 819 | ||||
| 820 | timebased_track->ddEdx_FDC = locdEdx_FDC; | |||
| 821 | timebased_track->ddx_FDC = locdx_FDC; | |||
| 822 | timebased_track->dNumHitsUsedFordEdx_FDC = locNumHitsUsedFordEdx_FDC; | |||
| 823 | timebased_track->ddEdx_CDC = locdEdx_CDC; | |||
| 824 | timebased_track->ddEdx_CDC_amp = locdEdx_CDC_amp; | |||
| 825 | timebased_track->ddx_CDC = locdx_CDC; | |||
| 826 | timebased_track->ddx_CDC_amp = locdx_CDC_amp; | |||
| 827 | timebased_track->dNumHitsUsedFordEdx_CDC = locNumHitsUsedFordEdx_CDC; | |||
| 828 | ||||
| 829 | timebased_track->potential_cdc_hits_on_track = fitter->GetNumPotentialCDCHits(); | |||
| 830 | timebased_track->potential_fdc_hits_on_track = fitter->GetNumPotentialFDCHits(); | |||
| 831 | ||||
| 832 | timebased_track->AddAssociatedObject(track); | |||
| 833 | _data.push_back(timebased_track); | |||
| 834 | ||||
| 835 | return true; | |||
| 836 | break; | |||
| 837 | } | |||
| 838 | case DTrackFitter::kFitSuccess: | |||
| 839 | { | |||
| 840 | // Create a new time-based track object | |||
| 841 | DTrackTimeBased *timebased_track = new DTrackTimeBased(); | |||
| 842 | *static_cast<DTrackingData*>(timebased_track) = fitter->GetFitParameters(); | |||
| 843 | ||||
| 844 | timebased_track->setTime(mStartTime); | |||
| 845 | timebased_track->chisq = fitter->GetChisq(); | |||
| 846 | timebased_track->Ndof = fitter->GetNdof(); | |||
| 847 | timebased_track->pulls = std::move(fitter->GetPulls()); | |||
| 848 | timebased_track->extrapolations=std::move(fitter->GetExtrapolations()); | |||
| 849 | timebased_track->IsSmoothed = fitter->GetIsSmoothed(); | |||
| 850 | timebased_track->trackid = track->id; | |||
| 851 | timebased_track->candidateid=track->candidateid; | |||
| 852 | timebased_track->flags=DTrackTimeBased::FLAG__GOODFIT; | |||
| 853 | ||||
| 854 | // Set the start time and add the list of start times | |||
| 855 | timebased_track->setT0(mStartTime,start_times[0].t0_sigma, mStartDetector); | |||
| 856 | timebased_track->start_times.assign(start_times.begin(), start_times.end()); | |||
| 857 | ||||
| 858 | if (DEBUG_HISTS){ | |||
| 859 | int id=0; | |||
| 860 | if (mStartDetector==SYS_CDC) id=1; | |||
| 861 | else if (mStartDetector==SYS_FDC) id=2; | |||
| 862 | else if (mStartDetector==SYS_BCAL) id=3; | |||
| 863 | else if (mStartDetector==SYS_FCAL) id=4; | |||
| 864 | else if (mStartDetector==SYS_TOF) id=5; | |||
| 865 | ||||
| 866 | Hstart_time->Fill(start_times[0].t0,id); | |||
| 867 | } | |||
| 868 | ||||
| 869 | ||||
| 870 | // Add hits used as associated objects | |||
| 871 | const vector<const DCDCTrackHit*> &cdchits = fitter->GetCDCFitHits(); | |||
| 872 | const vector<const DFDCPseudo*> &fdchits = fitter->GetFDCFitHits(); | |||
| 873 | ||||
| 874 | unsigned int num_fdc_potential=fitter->GetNumPotentialFDCHits(); | |||
| 875 | unsigned int num_cdc_potential=fitter->GetNumPotentialCDCHits(); | |||
| 876 | ||||
| 877 | DTrackTimeBased::hit_usage_t temp; | |||
| 878 | temp.inner_layer=0; | |||
| 879 | temp.outer_layer=0; | |||
| 880 | temp.total_hits=num_cdc_potential; | |||
| 881 | if (cdchits.size()>0){ | |||
| 882 | temp.inner_layer=cdchits[0]->wire->ring; | |||
| 883 | temp.outer_layer=cdchits[cdchits.size()-1]->wire->ring; | |||
| 884 | } | |||
| 885 | timebased_track->cdc_hit_usage=temp; | |||
| 886 | ||||
| 887 | // Reset the structure | |||
| 888 | temp.inner_layer=0; | |||
| 889 | temp.outer_layer=0; | |||
| 890 | temp.total_hits=num_fdc_potential; | |||
| 891 | if (fdchits.size()>0){ | |||
| 892 | temp.inner_layer=fdchits[0]->wire->layer; | |||
| 893 | temp.outer_layer=fdchits[fdchits.size()-1]->wire->layer; | |||
| 894 | } | |||
| 895 | timebased_track->fdc_hit_usage=temp; | |||
| 896 | ||||
| 897 | for(unsigned int m=0; m<cdchits.size(); m++) | |||
| 898 | timebased_track->AddAssociatedObject(cdchits[m]); | |||
| 899 | for(unsigned int m=0; m<fdchits.size(); m++) | |||
| 900 | timebased_track->AddAssociatedObject(fdchits[m]); | |||
| 901 | ||||
| 902 | timebased_track->measured_cdc_hits_on_track = cdchits.size(); | |||
| 903 | timebased_track->measured_fdc_hits_on_track = fdchits.size(); | |||
| 904 | ||||
| 905 | // dEdx | |||
| 906 | double locdEdx_FDC, locdx_FDC, locdEdx_CDC, locdEdx_CDC_amp; | |||
| 907 | double locdx_CDC,locdx_CDC_amp; | |||
| 908 | unsigned int locNumHitsUsedFordEdx_FDC, locNumHitsUsedFordEdx_CDC; | |||
| 909 | pid_algorithm->CalcDCdEdx(timebased_track, locdEdx_FDC, locdx_FDC, locdEdx_CDC, locdEdx_CDC_amp,locdx_CDC,locdx_CDC_amp, locNumHitsUsedFordEdx_FDC, locNumHitsUsedFordEdx_CDC); | |||
| 910 | ||||
| 911 | timebased_track->ddEdx_FDC = locdEdx_FDC; | |||
| 912 | timebased_track->ddx_FDC = locdx_FDC; | |||
| 913 | timebased_track->dNumHitsUsedFordEdx_FDC = locNumHitsUsedFordEdx_FDC; | |||
| 914 | timebased_track->ddEdx_CDC = locdEdx_CDC; | |||
| 915 | timebased_track->ddEdx_CDC_amp= locdEdx_CDC_amp; | |||
| 916 | timebased_track->ddx_CDC = locdx_CDC; | |||
| 917 | timebased_track->ddx_CDC_amp= locdx_CDC_amp; | |||
| 918 | timebased_track->dNumHitsUsedFordEdx_CDC = locNumHitsUsedFordEdx_CDC; | |||
| 919 | ||||
| 920 | // Set CDC ring & FDC plane hit patterns before candidate and wirebased tracks are associated | |||
| 921 | vector<const DCDCTrackHit*> tempCDCTrackHits; | |||
| 922 | vector<const DFDCPseudo*> tempFDCPseudos; | |||
| 923 | timebased_track->Get(tempCDCTrackHits); | |||
| 924 | timebased_track->Get(tempFDCPseudos); | |||
| 925 | timebased_track->dCDCRings = pid_algorithm->Get_CDCRingBitPattern(tempCDCTrackHits); | |||
| 926 | timebased_track->dFDCPlanes = pid_algorithm->Get_FDCPlaneBitPattern(tempFDCPseudos); | |||
| 927 | ||||
| 928 | timebased_track->potential_cdc_hits_on_track = fitter->GetNumPotentialCDCHits(); | |||
| 929 | timebased_track->potential_fdc_hits_on_track = fitter->GetNumPotentialFDCHits(); | |||
| 930 | ||||
| 931 | // Add DTrack object as associate object | |||
| 932 | timebased_track->AddAssociatedObject(track); | |||
| 933 | ||||
| 934 | // Compute the figure-of-merit based on tracking | |||
| 935 | timebased_track->FOM = TMath::Prob(timebased_track->chisq, timebased_track->Ndof); | |||
| 936 | //_DBG_<< "FOM: " << timebased_track->FOM << endl; | |||
| 937 | ||||
| 938 | _data.push_back(timebased_track); | |||
| 939 | ||||
| 940 | return true; | |||
| 941 | break; | |||
| 942 | ||||
| 943 | } | |||
| 944 | default: | |||
| 945 | break; | |||
| 946 | } | |||
| 947 | return false; | |||
| 948 | } | |||
| 949 | ||||
| 950 | ||||
| 951 | // Create a track with a mass hypothesis that was not present in the list of | |||
| 952 | // fitted tracks from an existing fitted track. | |||
| 953 | void DTrackTimeBased_factory::AddMissingTrackHypothesis(vector<DTrackTimeBased*>&tracks_to_add, | |||
| 954 | const DTrackTimeBased *src_track, | |||
| 955 | double my_mass, | |||
| 956 | double q, | |||
| 957 | JEventLoop *loop){ | |||
| 958 | ||||
| 959 | // Create a new time-based track object | |||
| 960 | DTrackTimeBased *timebased_track = new DTrackTimeBased(); | |||
| 961 | *static_cast<DTrackingData*>(timebased_track) = *static_cast<const DTrackingData*>(src_track); | |||
| 962 | ||||
| 963 | // Get the hits used in the fit | |||
| 964 | vector<const DCDCTrackHit *>src_cdchits; | |||
| 965 | src_track->GetT(src_cdchits); | |||
| 966 | vector<const DFDCPseudo *>src_fdchits; | |||
| 967 | src_track->GetT(src_fdchits); | |||
| 968 | ||||
| 969 | // Copy over DKinematicData part from the result of a successful fit | |||
| 970 | timebased_track->setPID(IDTrack(q, my_mass)); | |||
| 971 | timebased_track->chisq = src_track->chisq; | |||
| 972 | timebased_track->Ndof = src_track->Ndof; | |||
| 973 | timebased_track->pulls = src_track->pulls; | |||
| 974 | timebased_track->extrapolations = src_track->extrapolations; | |||
| 975 | timebased_track->trackid = src_track->id; | |||
| 976 | timebased_track->candidateid=src_track->candidateid; | |||
| 977 | timebased_track->FOM=src_track->FOM; | |||
| 978 | timebased_track->cdc_hit_usage=src_track->cdc_hit_usage; | |||
| 979 | timebased_track->fdc_hit_usage=src_track->fdc_hit_usage; | |||
| 980 | timebased_track->flags=DTrackTimeBased::FLAG__USED_OTHER_HYPOTHESIS; | |||
| 981 | ||||
| 982 | // Add list of start times | |||
| 983 | timebased_track->start_times.assign(src_track->start_times.begin(), | |||
| 984 | src_track->start_times.end()); | |||
| 985 | // Set the start time we used | |||
| 986 | timebased_track->setT0(timebased_track->start_times[0].t0, | |||
| 987 | timebased_track->start_times[0].t0_sigma, | |||
| 988 | timebased_track->start_times[0].system); | |||
| 989 | ||||
| 990 | // Add DTrack object as associate object | |||
| 991 | vector<const DTrackWireBased*>wire_based_track; | |||
| 992 | src_track->GetT(wire_based_track); | |||
| 993 | timebased_track->AddAssociatedObject(wire_based_track[0]); | |||
| 994 | ||||
| 995 | // (Partially) compensate for the difference in energy loss between the | |||
| 996 | // source track and a particle of mass my_mass | |||
| 997 | DVector3 position,momentum; | |||
| 998 | if (timebased_track->extrapolations.at(SYS_CDC).size()>0){ | |||
| 999 | unsigned int index=timebased_track->extrapolations.at(SYS_CDC).size()-1; | |||
| 1000 | position=timebased_track->extrapolations[SYS_CDC][index].position; | |||
| 1001 | momentum=timebased_track->extrapolations[SYS_CDC][index].momentum; | |||
| 1002 | } | |||
| 1003 | else if (timebased_track->extrapolations.at(SYS_FDC).size()>0){ | |||
| 1004 | unsigned int index=timebased_track->extrapolations.at(SYS_FDC).size()-1; | |||
| 1005 | position=timebased_track->extrapolations[SYS_FDC][index].position; | |||
| 1006 | momentum=timebased_track->extrapolations[SYS_FDC][index].momentum; | |||
| 1007 | } | |||
| 1008 | ||||
| 1009 | DTrackFitter::fit_status_t status=DTrackFitter::kFitNotDone; | |||
| 1010 | if (momentum.Mag()>0.){ | |||
| 1011 | CorrectForELoss(position,momentum,q,my_mass); | |||
| 1012 | timebased_track->setMomentum(momentum); | |||
| 1013 | timebased_track->setPosition(position); | |||
| 1014 | // Redo the fit with the new position and momentum as initial guesses | |||
| 1015 | fitter->Reset(); | |||
| 1016 | fitter->SetFitType(DTrackFitter::kTimeBased); | |||
| 1017 | status = fitter->FindHitsAndFitTrack(*timebased_track, | |||
| 1018 | timebased_track->extrapolations,loop, | |||
| 1019 | my_mass, | |||
| 1020 | src_cdchits.size()+2*src_fdchits.size(), | |||
| 1021 | timebased_track->t0(), | |||
| 1022 | timebased_track->t0_detector()); | |||
| 1023 | // if the fit returns chisq=-1, something went terribly wrong. Do not | |||
| 1024 | // update the parameters for the track... | |||
| 1025 | if (fitter->GetChisq()<0) status=DTrackFitter::kFitFailed; | |||
| 1026 | ||||
| 1027 | // if the fit flips the charge of the track, then this is bad as well | |||
| 1028 | if(q != fitter->GetFitParameters().charge()) | |||
| 1029 | status=DTrackFitter::kFitFailed; | |||
| 1030 | ||||
| 1031 | // if we can't refit the track, it is likely of poor quality, so stop here and do not add the hypothesis | |||
| 1032 | if(status == DTrackFitter::kFitFailed) { | |||
| 1033 | delete timebased_track; | |||
| 1034 | return; | |||
| 1035 | } | |||
| 1036 | ||||
| 1037 | if (status==DTrackFitter::kFitSuccess){ | |||
| 1038 | timebased_track->chisq = fitter->GetChisq(); | |||
| 1039 | timebased_track->Ndof = fitter->GetNdof(); | |||
| 1040 | timebased_track->pulls = std::move(fitter->GetPulls()); | |||
| 1041 | timebased_track->extrapolations=std::move(fitter->GetExtrapolations()); | |||
| 1042 | timebased_track->IsSmoothed = fitter->GetIsSmoothed(); | |||
| 1043 | *static_cast<DTrackingData*>(timebased_track) = fitter->GetFitParameters(); | |||
| 1044 | timebased_track->flags=DTrackTimeBased::FLAG__GOODFIT; | |||
| 1045 | ||||
| 1046 | timebased_track->setTime(timebased_track->start_times[0].t0); | |||
| 1047 | ||||
| 1048 | // Add hits used as associated objects | |||
| 1049 | const vector<const DCDCTrackHit*> &cdchits = fitter->GetCDCFitHits(); | |||
| 1050 | const vector<const DFDCPseudo*> &fdchits = fitter->GetFDCFitHits(); | |||
| 1051 | ||||
| 1052 | unsigned int num_fdc_potential=fitter->GetNumPotentialFDCHits(); | |||
| 1053 | unsigned int num_cdc_potential=fitter->GetNumPotentialCDCHits(); | |||
| 1054 | ||||
| 1055 | DTrackTimeBased::hit_usage_t temp; | |||
| 1056 | temp.inner_layer=0; | |||
| 1057 | temp.outer_layer=0; | |||
| 1058 | temp.total_hits=num_cdc_potential; | |||
| 1059 | if (cdchits.size()>0){ | |||
| 1060 | temp.inner_layer=cdchits[0]->wire->ring; | |||
| 1061 | temp.outer_layer=cdchits[cdchits.size()-1]->wire->ring; | |||
| 1062 | } | |||
| 1063 | timebased_track->cdc_hit_usage=temp; | |||
| 1064 | ||||
| 1065 | // Reset the structure | |||
| 1066 | temp.inner_layer=0; | |||
| 1067 | temp.outer_layer=0; | |||
| 1068 | temp.total_hits=num_fdc_potential; | |||
| 1069 | if (fdchits.size()>0){ | |||
| 1070 | temp.inner_layer=fdchits[0]->wire->layer; | |||
| 1071 | temp.outer_layer=fdchits[fdchits.size()-1]->wire->layer; | |||
| 1072 | } | |||
| 1073 | timebased_track->fdc_hit_usage=temp; | |||
| 1074 | ||||
| 1075 | for(unsigned int m=0; m<cdchits.size(); m++) | |||
| 1076 | timebased_track->AddAssociatedObject(cdchits[m]); | |||
| 1077 | for(unsigned int m=0; m<fdchits.size(); m++) | |||
| 1078 | timebased_track->AddAssociatedObject(fdchits[m]); | |||
| 1079 | ||||
| 1080 | timebased_track->measured_cdc_hits_on_track = cdchits.size(); | |||
| 1081 | timebased_track->measured_fdc_hits_on_track = fdchits.size(); | |||
| 1082 | ||||
| 1083 | // Compute the figure-of-merit based on tracking | |||
| 1084 | timebased_track->FOM = TMath::Prob(timebased_track->chisq, timebased_track->Ndof); | |||
| 1085 | ||||
| 1086 | } | |||
| 1087 | } | |||
| 1088 | ||||
| 1089 | if (status!=DTrackFitter::kFitSuccess){ | |||
| 1090 | for(unsigned int m=0; m<src_fdchits.size(); m++) | |||
| 1091 | timebased_track->AddAssociatedObject(src_fdchits[m]); | |||
| 1092 | for(unsigned int m=0; m<src_cdchits.size(); m++) | |||
| 1093 | timebased_track->AddAssociatedObject(src_cdchits[m]); | |||
| 1094 | ||||
| 1095 | timebased_track->measured_cdc_hits_on_track = src_cdchits.size(); | |||
| 1096 | timebased_track->measured_fdc_hits_on_track = src_fdchits.size(); | |||
| 1097 | } | |||
| 1098 | ||||
| 1099 | // dEdx | |||
| 1100 | double locdEdx_FDC, locdx_FDC, locdEdx_CDC, locdEdx_CDC_amp; | |||
| 1101 | double locdx_CDC,locdx_CDC_amp; | |||
| 1102 | unsigned int locNumHitsUsedFordEdx_FDC, locNumHitsUsedFordEdx_CDC; | |||
| 1103 | pid_algorithm->CalcDCdEdx(timebased_track, locdEdx_FDC, locdx_FDC, locdEdx_CDC, locdEdx_CDC_amp,locdx_CDC,locdx_CDC_amp, locNumHitsUsedFordEdx_FDC, locNumHitsUsedFordEdx_CDC); | |||
| 1104 | ||||
| 1105 | timebased_track->ddEdx_FDC = locdEdx_FDC; | |||
| 1106 | timebased_track->ddx_FDC = locdx_FDC; | |||
| 1107 | timebased_track->dNumHitsUsedFordEdx_FDC = locNumHitsUsedFordEdx_FDC; | |||
| 1108 | timebased_track->ddEdx_CDC = locdEdx_CDC; | |||
| 1109 | timebased_track->ddEdx_CDC_amp = locdEdx_CDC_amp; | |||
| 1110 | timebased_track->ddx_CDC = locdx_CDC; | |||
| 1111 | timebased_track->ddx_CDC_amp = locdx_CDC_amp; | |||
| 1112 | timebased_track->dNumHitsUsedFordEdx_CDC = locNumHitsUsedFordEdx_CDC; | |||
| 1113 | ||||
| 1114 | // Set CDC ring & FDC plane hit patterns before candidate and wirebased tracks are associated | |||
| 1115 | vector<const DCDCTrackHit*> tempCDCTrackHits; | |||
| 1116 | vector<const DFDCPseudo*> tempFDCPseudos; | |||
| 1117 | timebased_track->Get(tempCDCTrackHits); | |||
| 1118 | timebased_track->Get(tempFDCPseudos); | |||
| 1119 | timebased_track->dCDCRings = pid_algorithm->Get_CDCRingBitPattern(tempCDCTrackHits); | |||
| 1120 | timebased_track->dFDCPlanes = pid_algorithm->Get_FDCPlaneBitPattern(tempFDCPseudos); | |||
| 1121 | ||||
| 1122 | timebased_track->potential_cdc_hits_on_track = fitter->GetNumPotentialCDCHits(); | |||
| 1123 | timebased_track->potential_fdc_hits_on_track = fitter->GetNumPotentialFDCHits(); | |||
| 1124 | ||||
| 1125 | tracks_to_add.push_back(timebased_track); | |||
| 1126 | } | |||
| 1127 | ||||
| 1128 | ||||
| 1129 | // If the fit failed for certain hypotheses, fill in the gaps using data from | |||
| 1130 | // successful fits for each candidate. | |||
| 1131 | bool DTrackTimeBased_factory::InsertMissingHypotheses(JEventLoop *loop){ | |||
| 1132 | if (_data.size()==0) return false; | |||
| 1133 | ||||
| 1134 | // Make sure the tracks are ordered by candidate id | |||
| 1135 | sort(_data.begin(),_data.end(),DTrackTimeBased_cmp); | |||
| 1136 | ||||
| 1137 | JObject::oid_t old_id=_data[0]->candidateid; | |||
| 1138 | unsigned int mass_bits=0; | |||
| 1139 | double q=_data[0]->charge(); | |||
| 1140 | bool flipped_charge=false; | |||
| 1141 | vector<DTrackTimeBased*>myhypotheses; | |||
| 1142 | vector<DTrackTimeBased*>tracks_to_add; | |||
| 1143 | for (size_t i=0;i<_data.size();i++){ | |||
| 1144 | if (_data[i]->candidateid!=old_id){ | |||
| 1145 | AddMissingTrackHypotheses(mass_bits,tracks_to_add,myhypotheses,q, | |||
| 1146 | flipped_charge,loop); | |||
| 1147 | ||||
| 1148 | // Clear the myhypotheses vector for the next track | |||
| 1149 | myhypotheses.clear(); | |||
| 1150 | // Reset flags and charge | |||
| 1151 | q=_data[i]->charge(); | |||
| 1152 | flipped_charge=false; | |||
| 1153 | // Set the bit for this mass hypothesis | |||
| 1154 | mass_bits = 1<<_data[i]->PID(); | |||
| 1155 | ||||
| 1156 | // Add the data to the myhypotheses vector | |||
| 1157 | myhypotheses.push_back(_data[i]); | |||
| 1158 | } | |||
| 1159 | else{ | |||
| 1160 | myhypotheses.push_back(_data[i]); | |||
| 1161 | ||||
| 1162 | // Set the bit for this mass hypothesis | |||
| 1163 | mass_bits |= 1<< _data[i]->PID(); | |||
| 1164 | ||||
| 1165 | // Check if the sign of the charge has flipped | |||
| 1166 | if (_data[i]->charge()!=q) flipped_charge=true; | |||
| 1167 | } | |||
| 1168 | ||||
| 1169 | old_id=_data[i]->candidateid; | |||
| 1170 | } | |||
| 1171 | // Deal with last track candidate | |||
| 1172 | AddMissingTrackHypotheses(mass_bits,tracks_to_add,myhypotheses,q, | |||
| 1173 | flipped_charge,loop); | |||
| 1174 | ||||
| 1175 | // Add the new list of tracks to the output list | |||
| 1176 | if (tracks_to_add.size()>0){ | |||
| 1177 | for (size_t i=0;i<tracks_to_add.size();i++){ | |||
| 1178 | _data.push_back(tracks_to_add[i]); | |||
| 1179 | } | |||
| 1180 | // Make sure the tracks are ordered by candidate id | |||
| 1181 | sort(_data.begin(),_data.end(),DTrackTimeBased_cmp); | |||
| 1182 | } | |||
| 1183 | ||||
| 1184 | return true; | |||
| 1185 | } | |||
| 1186 | ||||
| 1187 | // Use the FastSwim method in DReferenceTrajectory to propagate back to the | |||
| 1188 | // POCA to the beam line, adding a bit of energy at each step that would have | |||
| 1189 | // been lost had the particle emerged from the target. | |||
| 1190 | void DTrackTimeBased_factory::CorrectForELoss(DVector3 &position,DVector3 &momentum,double q,double my_mass){ | |||
| 1191 | DReferenceTrajectory::swim_step_t dummy_step; | |||
| 1192 | DReferenceTrajectory rt(fitter->GetDMagneticFieldMap(),q,&dummy_step); | |||
| 1193 | rt.SetDGeometry(geom); | |||
| 1194 | rt.SetMass(my_mass); | |||
| 1195 | rt.SetPLossDirection(DReferenceTrajectory::kBackward); | |||
| 1196 | DVector3 last_pos,last_mom; | |||
| 1197 | DVector3 origin(0.,0.,65.); | |||
| 1198 | DVector3 dir(0.,0.,1.); | |||
| 1199 | rt.FastSwim(position,momentum,last_pos,last_mom,q,origin,dir,300.); | |||
| 1200 | position=last_pos; | |||
| 1201 | momentum=last_mom; | |||
| 1202 | } | |||
| 1203 | ||||
| 1204 | // Fill in all missing hypotheses for a given track candidate | |||
| 1205 | void DTrackTimeBased_factory::AddMissingTrackHypotheses(unsigned int mass_bits, | |||
| 1206 | vector<DTrackTimeBased*>&tracks_to_add, | |||
| 1207 | vector<DTrackTimeBased *>&myhypotheses, | |||
| 1208 | double q, | |||
| 1209 | bool flipped_charge, | |||
| 1210 | JEventLoop *loop){ | |||
| 1211 | ||||
| 1212 | unsigned int last_index=myhypotheses.size()-1; | |||
| 1213 | unsigned int index=0; | |||
| 1214 | if (q>0){ | |||
| 1215 | if (flipped_charge){ | |||
| 1216 | /*if ((mass_bits & (1<<AntiProton))==0){ | |||
| 1217 | AddMissingTrackHypothesis(tracks_to_add,myhypotheses[last_index], | |||
| 1218 | ParticleMass(Proton),-1.,loop); | |||
| 1219 | } | |||
| 1220 | */ | |||
| 1221 | for (unsigned int i=0;i<mass_hypotheses_negative.size();i++){ | |||
| 1222 | if ((mass_bits & (1<<mass_hypotheses_negative[i]))==0){ | |||
| 1223 | if (mass_hypotheses_negative[i]>13) index=last_index; | |||
| 1224 | else index=0; | |||
| 1225 | AddMissingTrackHypothesis(tracks_to_add,myhypotheses[index], | |||
| 1226 | ParticleMass(Particle_t(mass_hypotheses_negative[i])), | |||
| 1227 | -1.,loop); | |||
| 1228 | } | |||
| 1229 | } | |||
| 1230 | } | |||
| 1231 | /* if ((mass_bits & (1<<Proton))==0){ | |||
| 1232 | AddMissingTrackHypothesis(tracks_to_add,myhypotheses[last_index], | |||
| 1233 | ParticleMass(Proton),+1.,loop); | |||
| 1234 | } */ | |||
| 1235 | for (unsigned int i=0;i<mass_hypotheses_positive.size();i++){ | |||
| 1236 | if ((mass_bits & (1<<mass_hypotheses_positive[i]))==0){ | |||
| 1237 | if (mass_hypotheses_positive[i]>13) index=last_index; | |||
| 1238 | else index=0; | |||
| 1239 | AddMissingTrackHypothesis(tracks_to_add,myhypotheses[index], | |||
| 1240 | ParticleMass(Particle_t(mass_hypotheses_positive[i])), | |||
| 1241 | +1.,loop); | |||
| 1242 | } | |||
| 1243 | } | |||
| 1244 | } | |||
| 1245 | else{ | |||
| 1246 | /* | |||
| 1247 | if ((mass_bits & (1<<AntiProton))==0){ | |||
| 1248 | AddMissingTrackHypothesis(tracks_to_add,myhypotheses[last_index], | |||
| 1249 | ParticleMass(Proton),-1.,loop); | |||
| 1250 | } */ | |||
| 1251 | for (unsigned int i=0;i<mass_hypotheses_negative.size();i++){ | |||
| 1252 | if ((mass_bits & (1<<mass_hypotheses_negative[i]))==0){ | |||
| 1253 | if (mass_hypotheses_negative[i]>13) index=last_index; | |||
| 1254 | else index=0; | |||
| 1255 | AddMissingTrackHypothesis(tracks_to_add,myhypotheses[index], | |||
| 1256 | ParticleMass(Particle_t(mass_hypotheses_negative[i])), | |||
| 1257 | -1.,loop); | |||
| 1258 | } | |||
| 1259 | } | |||
| 1260 | if (flipped_charge){ | |||
| 1261 | /*if ((mass_bits & (1<<Proton))==0){ | |||
| 1262 | AddMissingTrackHypothesis(tracks_to_add,myhypotheses[last_index], | |||
| 1263 | ParticleMass(Proton),+1.,loop); | |||
| 1264 | } */ | |||
| 1265 | for (unsigned int i=0;i<mass_hypotheses_positive.size();i++){ | |||
| 1266 | if ((mass_bits & (1<<mass_hypotheses_positive[i]))==0){ | |||
| 1267 | if (mass_hypotheses_positive[i]>13) index=last_index; | |||
| 1268 | else index=0; | |||
| 1269 | AddMissingTrackHypothesis(tracks_to_add,myhypotheses[index], | |||
| 1270 | ParticleMass(Particle_t(mass_hypotheses_positive[i])), | |||
| 1271 | +1.,loop); | |||
| 1272 | } | |||
| 1273 | } | |||
| 1274 | } | |||
| 1275 | } | |||
| 1276 | } |
| 1 | // $Id: JObject.h 1709 2006-04-26 20:34:03Z davidl $ |
| 2 | // |
| 3 | // File: JObject.h |
| 4 | // Created: Wed Aug 17 10:57:09 EDT 2005 |
| 5 | // Creator: davidl (on Darwin wire129.jlab.org 7.8.0 powerpc) |
| 6 | // |
| 7 | |
| 8 | #ifndef _JObject_ |
| 9 | #define _JObject_ |
| 10 | |
| 11 | #include <cstdio> |
| 12 | #include <sstream> |
| 13 | #include <cassert> |
| 14 | #include <map> |
| 15 | #include <vector> |
| 16 | #include <set> |
| 17 | #include <string> |
| 18 | #include <cstdint> |
| 19 | #include <type_traits> |
| 20 | #include <typeinfo> |
| 21 | #include <stdint.h> |
| 22 | using std::pair; |
| 23 | using std::map; |
| 24 | using std::set; |
| 25 | using std::vector; |
| 26 | using std::string; |
| 27 | using std::stringstream; |
| 28 | |
| 29 | // The following is here just so we can use ROOT's THtml class to generate documentation. |
| 30 | #include "cint.h" |
| 31 | |
| 32 | |
| 33 | /// The JObject class is a base class for all data classes. |
| 34 | /// (See JFactory and JFactory_base for algorithm classes.) |
| 35 | /// |
| 36 | /// |
| 37 | /// The following line should be included in the public definition of all |
| 38 | /// classes which inherit from JObject with the argument being the name of the |
| 39 | /// class without any quotes. |
| 40 | /// e.g. for a class named "MyClass" : |
| 41 | /// |
| 42 | /// public: |
| 43 | /// JOBJECT_PUBLIC(MyClass); |
| 44 | /// |
| 45 | /// This will define a virtual method <i>className()</i> and a static |
| 46 | /// method <i>static_className()</i> that are used by JANA to identify |
| 47 | /// the object's last generation in the inheritance chain by name. This |
| 48 | /// also allows for possible upgrades to JANA in the future without |
| 49 | /// requiring classes that inherit from JObject to be redefined explicity. |
| 50 | #define JOBJECT_PUBLIC(T)virtual const char* className(void) const {return static_className ();} static const char* static_className(void) {return "T";} virtual JObject* Clone() const {return CloneObject<T>( *this ) ;} \ |
| 51 | virtual const char* className(void) const {return static_className();} \ |
| 52 | static const char* static_className(void) {return #T;} \ |
| 53 | virtual JObject* Clone() const {return CloneObject<T>( *this );} |
| 54 | |
| 55 | |
| 56 | |
| 57 | // Place everything in JANA namespace |
| 58 | namespace jana{ |
| 59 | |
| 60 | class JFactory_base; |
| 61 | class JEventLoop; |
| 62 | |
| 63 | class JObject{ |
| 64 | |
| 65 | public: |
| 66 | JOBJECT_PUBLIC(JObject)virtual const char* className(void) const {return static_className ();} static const char* static_className(void) {return "JObject" ;} virtual JObject* Clone() const {return CloneObject<JObject >( *this );}; |
| 67 | |
| 68 | typedef unsigned long long oid_t; |
| 69 | |
| 70 | JObject() : id((oid_t)this),append_types(false),factory(NULL__null) {} |
| 71 | JObject( oid_t aId ) : id( aId ),append_types(false),factory(NULL__null) {} |
| 72 | |
| 73 | virtual ~JObject(){ |
| 74 | for(unsigned int i=0; i<auto_delete.size(); i++)delete auto_delete[i]; |
| 75 | } |
| 76 | |
| 77 | // Copy constructor |
| 78 | JObject(const JObject& o) : |
| 79 | id( (oid_t)this ), append_types(o.append_types), associated( |
| 80 | o.associated), auto_delete(), messagelog( |
| 81 | o.messagelog), factory(o.factory) { |
| 82 | // Deep copy the objects in auto_delete vector |
| 83 | for( auto obj : o.auto_delete ) { |
| 84 | auto_delete.push_back( obj->Clone() ); |
| 85 | } |
| 86 | } |
| 87 | |
| 88 | // Move constructor |
| 89 | JObject( const JObject&& o ) : id( (oid_t)this ), append_types(o.append_types), associated( |
| 90 | std::move(o.associated)), auto_delete(std::move(o.auto_delete)), messagelog( std::move( |
| 91 | o.messagelog)), factory(o.factory) { |
| 92 | factory = nullptr; |
| 93 | } |
| 94 | |
| 95 | // Assignment operator |
| 96 | JObject& operator=( const JObject& o) { |
| 97 | if( this == &o ) return *this; |
| 98 | append_types = o.append_types; |
| 99 | associated = o.associated; |
| 100 | messagelog = o.messagelog; |
| 101 | factory = o.factory; |
| 102 | auto_delete.clear(); |
| 103 | for( auto obj : o.auto_delete ) { |
| 104 | auto_delete.push_back( obj->Clone() ); |
| 105 | } |
| 106 | return *this; |
| 107 | } |
| 108 | |
| 109 | // Define template static method for cloning the object of this type. This would get called from a |
| 110 | // virtual method define for each subclass to polymorphically clone the objects. |
| 111 | template<typename TYPE> |
| 112 | static typename std::enable_if< |
| 113 | (std::is_abstract <TYPE>::value || !std::is_copy_constructible<TYPE>::value), JObject*>::type CloneObject( |
| 114 | const TYPE& obj) { |
| 115 | // For abstract method it will return null pointer, and probably will cause problems. But |
| 116 | // there sohuld not be any object of abstract class to start with. |
| 117 | return nullptr; |
| 118 | } |
| 119 | template<typename TYPE> |
| 120 | static typename std::enable_if< |
| 121 | (!std::is_abstract <TYPE>::value |
| 122 | && std::is_copy_constructible <TYPE>::value ), JObject*>::type CloneObject( |
| 123 | const TYPE& obj) { |
| 124 | // For non-abstract class create a new object using the copy constructor of the class |
| 125 | // given in the template argument. |
| 126 | return new TYPE(obj); |
| 127 | } |
| 128 | |
| 129 | /// Test if this object is of type T by checking its className() against T::static_className() |
| 130 | template<typename T> bool IsA(const T *t) const {return dynamic_cast<const T*>(this)!=0L;} |
| 131 | |
| 132 | /// Test if this object is of type T (or a descendant) by attempting a dynamic cast |
| 133 | template<typename T> bool IsAT(const T *t) const {return dynamic_cast<const T*>(this)!=0L;} |
| 134 | |
| 135 | // Methods for handling associated objects |
| 136 | inline void AddAssociatedObject(const JObject *obj); |
| 137 | inline void AddAssociatedObjectAutoDelete(JObject *obj, bool auto_delete=true); |
| 138 | inline void RemoveAssociatedObject(const JObject *obj); |
| 139 | inline void ClearAssociatedObjects(void); |
| 140 | inline bool IsAssociated(const JObject* locObject) const {return (associated.find(locObject) != associated.end());} |
| 141 | template<typename T> void Get(vector<const T*> &ptrs, string classname="", int max_depth=1000000) const ; |
| 142 | template<typename T> void GetT(vector<const T*> &ptrs) const ; |
| 143 | template<typename T> void GetSingle(const T* &ptrs, string classname="") const ; |
| 144 | template<typename T> void GetSingleT(const T* &ptrs) const ; |
| 145 | template<typename T> void GetAssociatedAncestors(set<const JObject*> &already_checked, int &max_depth, set<const T*> &objs_found, string classname="") const; |
| 146 | template<typename T> void GetAssociatedDescendants(JEventLoop *loop, vector<const T*> &associatedTo, int max_depth=1000000); |
| 147 | void GetAssociatedDescendants(JEventLoop *loop, vector<const JObject*> &associatedTo, int max_depth=1000000); |
| 148 | |
| 149 | template<typename T,typename S> void CopyToVector(T itbegin, T itend, vector<const S*> &v) const; |
| 150 | |
| 151 | // Methods for handling pretty formatting for dumping to the screen or file |
| 152 | virtual void toStrings(vector<pair<string,string> > &items)const; |
| 153 | template<typename T> void AddString(vector<pair<string,string> > &items, const char *name, const char *format, const T &val) const; |
| 154 | |
| 155 | // Methods for attaching and retrieving log messages to/from object |
| 156 | void AddLog(string &message) const {messagelog.push_back(message);} |
| 157 | void AddLog(vector<string> &messages) const {messagelog.insert(messagelog.end(), messages.begin(), messages.end());} |
| 158 | void GetLog(vector<string> &messagelog) const {messagelog = this->messagelog;} |
| 159 | |
| 160 | // Misc methods |
| 161 | bool GetAppendTypes(void) const {return append_types;} ///< Get state of append_types flag (for AddString) |
| 162 | void SetAppendTypes(bool append_types){this->append_types=append_types;} ///< Set state of append_types flag (for AddString) |
| 163 | void SetFactoryPointer(JFactory_base *factory){this->factory=factory;} |
| 164 | JFactory_base * GetFactoryPointer(void) const {return factory;} |
| 165 | string GetName(void) const {return string(className());} |
| 166 | string GetTag(void) const ; |
| 167 | string GetNameTag(void) const {return GetName() + (GetTag()=="" ? "":":") + GetTag();} |
| 168 | |
| 169 | oid_t id; |
| 170 | |
| 171 | private: |
| 172 | |
| 173 | bool append_types; |
| 174 | set<const JObject*> associated; |
| 175 | // map<const JObject*, string> associated; replaced with set in jana 0.7.7 |
| 176 | vector<JObject*> auto_delete; |
| 177 | mutable vector<string> messagelog; |
| 178 | JFactory_base *factory; |
| 179 | |
| 180 | }; |
| 181 | |
| 182 | #if !defined(__CINT__) && !defined(__CLING__) |
| 183 | |
| 184 | |
| 185 | //-------------------------- |
| 186 | // AddAssociatedObject |
| 187 | //-------------------------- |
| 188 | void JObject::AddAssociatedObject(const JObject *obj) |
| 189 | { |
| 190 | /// Add a JObject to the list of associated objects |
| 191 | |
| 192 | assert(obj!=NULL)((obj!=__null) ? static_cast<void> (0) : __assert_fail ( "obj!=__null", "/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JObject.h" , 192, __PRETTY_FUNCTION__)); |
| 193 | |
| 194 | associated.insert(obj); |
| 195 | //associated[obj] = obj->className(); |
| 196 | } |
| 197 | |
| 198 | //-------------------------- |
| 199 | // AddAssociatedObjectAutoDelete |
| 200 | //-------------------------- |
| 201 | void JObject::AddAssociatedObjectAutoDelete(JObject *obj, bool auto_delete) |
| 202 | { |
| 203 | /// Add a JObject to the list of associated objects. If the auto_delete |
| 204 | /// flag is true, then automatically delete it when this object is |
| 205 | /// deleted. Otherwise, this behaves identically to the AddAssociatedObject |
| 206 | /// method. |
| 207 | /// |
| 208 | /// Note that if the object is removed via RemoveAssociatedObject(...) |
| 209 | /// then the object is NOT deleted. BUT, if the entire list of associated |
| 210 | /// objects is cleared via ClearAssociatedObjects, then the object will |
| 211 | /// be deleted. |
| 212 | |
| 213 | AddAssociatedObject(obj); |
| 214 | |
| 215 | if(auto_delete)this->auto_delete.push_back(obj); |
| 216 | } |
| 217 | |
| 218 | //-------------------------- |
| 219 | // RemoveAssociatedObject |
| 220 | //-------------------------- |
| 221 | void JObject::RemoveAssociatedObject(const JObject *obj) |
| 222 | { |
| 223 | /// Remove the specified JObject from the list of associated |
| 224 | /// objects. This will NOT delete the object even if the |
| 225 | /// object was added with the AddAssociatedObjectAutoDelete(...) |
| 226 | /// method with the auto_delete flag set. |
| 227 | |
| 228 | // map<const JObject*, string>::iterator iter = associated.find(obj); |
| 229 | auto iter = associated.find(obj); |
| 230 | |
| 231 | if(iter!=associated.end()){ |
| 232 | associated.erase(iter); |
| 233 | } |
| 234 | } |
| 235 | |
| 236 | //-------------------------- |
| 237 | // ClearAssociatedObjects |
| 238 | //-------------------------- |
| 239 | void JObject::ClearAssociatedObjects(void) |
| 240 | { |
| 241 | /// Remove all associated objects from the associated objects list. |
| 242 | /// This will also delete any objects that were added via the |
| 243 | /// AddAssociatedObjectAutoDelete(...) method with the auto_delete |
| 244 | /// flag set. |
| 245 | |
| 246 | // Clear pointers to associated objects |
| 247 | associated.clear(); |
| 248 | |
| 249 | // Delete objects in the auto_delete list |
| 250 | for(unsigned int i=0; i<auto_delete.size(); i++)delete auto_delete[i]; |
| 251 | auto_delete.clear(); |
| 252 | } |
| 253 | |
| 254 | //-------------------------- |
| 255 | // CopyToVector |
| 256 | // |
| 257 | // This litte utility method is here because the g++ 4.4.7 |
| 258 | // compiler was complaining about the declaration of |
| 259 | // set<const T*>::iterator it; in the Get method. Very strange |
| 260 | // but this at least avoids ever having to declare the variable. |
| 261 | //-------------------------- |
| 262 | template<typename T,typename S> |
| 263 | void JObject::CopyToVector(T itbegin, T itend, vector<const S*> &v) const |
| 264 | { |
| 265 | for(T it=itbegin; it!=itend; it++) v.push_back(*it); |
| 266 | } |
| 267 | |
| 268 | //-------------------------- |
| 269 | // Get |
| 270 | //-------------------------- |
| 271 | template<typename T> |
| 272 | void JObject::Get(vector<const T*> &ptrs, string classname, int max_depth) const |
| 273 | { |
| 274 | /// Fill the given vector with pointers to the associated objects of the |
| 275 | /// type on which the vector is based. The objects are chosen by matching |
| 276 | /// their class names (obtained via JObject::className()) either to the |
| 277 | /// one provided in classname or to T::static_className() if classname is |
| 278 | /// an empty string. Associations will be searched to a level of max_depth |
| 279 | /// to find all objects of the requested type. By default, max_depth is |
| 280 | /// set to a very large number so that all associations are found. To |
| 281 | /// limit the search to only objects directly associated with this one, |
| 282 | /// set max_depth to either "0" or "1". |
| 283 | /// |
| 284 | /// The contents of ptrs are cleared upon entry. |
| 285 | |
| 286 | if(classname=="")classname=T::static_className(); |
| 287 | |
| 288 | // Use the GetAssociatedAncestors method which may call itself |
| 289 | // recursively to search all levels of association (at or |
| 290 | // below this object. Objects for which this is an associated |
| 291 | // object are not checked for). |
| 292 | set<const JObject*> already_checked; |
| 293 | set<const T*> objs_found; |
| 294 | int my_max_depth = max_depth; |
| 295 | GetAssociatedAncestors(already_checked, my_max_depth, objs_found, classname); |
| 296 | |
| 297 | // Copy results into caller's container |
| 298 | ptrs.clear(); |
| 299 | CopyToVector(objs_found.begin(), objs_found.end(), ptrs); |
| 300 | // set<const T*>::iterator it; |
| 301 | // for(it=objs_found.begin(); it!=objs_found.end(); it++){ |
| 302 | // ptrs.push_back(*it); |
| 303 | // } |
| 304 | } |
| 305 | |
| 306 | //-------------------------- |
| 307 | // GetAssociatedAncestors |
| 308 | //-------------------------- |
| 309 | template<typename T> |
| 310 | void JObject::GetAssociatedAncestors(set<const JObject*> &already_checked, int &max_depth, set<const T*> &objs_found, string classname) const |
| 311 | { |
| 312 | /// Get associated objects of the specified type (either "T" or classname). |
| 313 | /// Check also for associated objects of any associated objects |
| 314 | /// to a level of max_depth associations. This method calls itself |
| 315 | /// recursively so care is taken to only check the associated objects |
| 316 | /// of each object encountered only once. |
| 317 | /// |
| 318 | /// The "already_checked" parameter should be passed in as an empty container |
| 319 | /// that is used to keep track of which objects had their direct associations |
| 320 | /// checked. "max_depth" indicates the maximum level of associations to check |
| 321 | /// (n.b. both "0" and "1" means only check direct associations.) This must |
| 322 | /// be passed as a reference to an existing int since it is modified in order |
| 323 | /// to keep track of the current depth in the recursive calls. Set max_depth |
| 324 | /// to a very high number (like 1000000) to check all associations. The |
| 325 | /// "objs_found" container will contain the actual associated objects found. |
| 326 | /// The objects are chosen by matching their class names (obtained via |
| 327 | /// JObject::className()) either to the one provided in "classname" or to |
| 328 | /// T::static_className() if classname is an empty string. |
| 329 | |
| 330 | if(already_checked.find(this) == already_checked.end()) already_checked.insert(this); |
| 331 | |
| 332 | if(classname=="")classname=T::static_className(); |
| 333 | max_depth--; |
| 334 | |
| 335 | //map<const JObject*, string>::const_iterator iter = associated.begin(); |
| 336 | for( auto obj : associated ){ |
| 337 | |
| 338 | // Add to list if appropriate |
| 339 | if( classname == obj->className() ){ |
| 340 | objs_found.insert( dynamic_cast<const T*>(obj) ); |
| 341 | } |
| 342 | |
| 343 | // Check this object's associated objects if appropriate |
| 344 | if(max_depth<=0) continue; |
| 345 | if(already_checked.find(obj) != already_checked.end()) continue; |
| 346 | already_checked.insert(obj); |
| 347 | obj->GetAssociatedAncestors(already_checked, max_depth, objs_found, classname); |
| 348 | } |
| 349 | |
| 350 | max_depth++; |
| 351 | } |
| 352 | |
| 353 | //-------------------------- |
| 354 | // GetAssociatedDescendants |
| 355 | //-------------------------- |
| 356 | template<typename T> |
| 357 | void GetAssociatedDescendants(JEventLoop *loop, vector<const T*> &associatedTo, int max_depth=1000000) |
| 358 | { |
| 359 | /// Find objects of type "T" for which this object appears in its |
| 360 | /// associated ancestors list. (This is kind of the opposite of |
| 361 | /// the "Get()" method.) |
| 362 | /// |
| 363 | /// WARNING: this must build and search the ancestor list of EVERY |
| 364 | /// object produced by EVERY factory. It is an expensive method |
| 365 | /// to call. Use it with great caution! |
| 366 | /// |
| 367 | /// WARNING: this only searches objects that have already been |
| 368 | /// created. It will not activate factories they may eventually |
| 369 | /// claim this as an associated object so the list returned may |
| 370 | /// be incomplete. |
| 371 | /// |
| 372 | /// WARNING: this templated method works by first calling the |
| 373 | /// JObject form and then dynamically casting each of those to see |
| 374 | /// if they of type "T". This makes this an even more expensive |
| 375 | /// call. Again, use with great caution!! |
| 376 | |
| 377 | vector<const JObject*> ajobjs; |
| 378 | GetAssociatedDescendants(loop, ajobjs, max_depth); |
| 379 | for(uint32_t i=0; i<ajobjs.size(); i++){ |
| 380 | const T *ptr = dynamic_cast<const T*>(ajobjs[i]); |
| 381 | if(ptr != NULL__null) associatedTo.push_back(ptr); |
| 382 | } |
| 383 | } |
| 384 | |
| 385 | //-------------------------- |
| 386 | // GetT |
| 387 | //-------------------------- |
| 388 | template<typename T> |
| 389 | void JObject::GetT(vector<const T*> &ptrs) const |
| 390 | { |
| 391 | /// Fill the given vector with pointers to the associated |
| 392 | /// JObjects of the type on which the vector is based. This is |
| 393 | /// similar to the Get() method except objects are selected |
| 394 | /// by attempting a dynamic_cast to type const T*. This allows |
| 395 | /// one to select a list of all objects who have a type T |
| 396 | /// somewhere in their inheritance chain. |
| 397 | /// |
| 398 | /// A potential issue with this method is that the dynamic_cast |
| 399 | /// does not always work correctly for objects created via a |
| 400 | /// plugin when the cast occurs outside of the plugin or |
| 401 | /// vice versa. |
| 402 | /// |
| 403 | /// The contents of ptrs are cleared upon entry. |
| 404 | |
| 405 | ptrs.clear(); |
| 406 | |
| 407 | //map<const JObject*, string>::const_iterator iter = associated.begin(); |
| 408 | //for(; iter!=associated.end(); iter++){ |
| 409 | for( auto obj : associated ){ |
| 410 | const T *ptr = dynamic_cast<const T*>(obj); |
| 411 | if(ptr != NULL__null)ptrs.push_back(ptr); |
| 412 | } |
| 413 | } |
| 414 | |
| 415 | //------------- |
| 416 | // GetSingle |
| 417 | //------------- |
| 418 | template<class T> |
| 419 | void JObject::GetSingle(const T* &t, string classname) const |
| 420 | { |
| 421 | /// This is a convenience method that can be used to get a pointer to the single |
| 422 | /// associate object of type T. |
| 423 | /// |
| 424 | /// The objects are chosen by matching their class names |
| 425 | /// (obtained via JObject::className()) either |
| 426 | /// to the one provided in classname or to T::static_className() |
| 427 | /// if classname is an empty string. |
| 428 | /// |
| 429 | /// If no object of the specified type is found, a NULL pointer is |
| 430 | /// returned. |
| 431 | |
| 432 | t = NULL__null; |
| 433 | |
| 434 | if(classname=="")classname=T::static_className(); |
| 435 | |
| 436 | //map<const JObject*, string>::const_iterator iter = associated.begin(); |
| 437 | //for(; iter!=associated.end(); iter++){ |
| 438 | for( auto obj : associated ){ |
| 439 | if( classname == obj->className() ){ |
| 440 | t = dynamic_cast<const T*>(obj); |
| 441 | if(t!=NULL__null)return; |
| 442 | } |
| 443 | } |
| 444 | } |
| 445 | |
| 446 | //------------- |
| 447 | // GetSingleT |
| 448 | //------------- |
| 449 | template<class T> |
| 450 | void JObject::GetSingleT(const T* &t) const |
| 451 | { |
| 452 | /// This is a convenience method that can be used to get a pointer to the single |
| 453 | /// associate object of type T. |
| 454 | /// |
| 455 | /// This is similar to the GetSingle() method except objects are selected |
| 456 | /// by attempting a dynamic_cast to type const T*. This allows |
| 457 | /// one to select a list of all objects who have a type T |
| 458 | /// somewhere in their inheritance chain. |
| 459 | /// The objects are chosen by matching their class names |
| 460 | /// (obtained via JObject::className()) either |
| 461 | /// to the one provided in classname or to T::static_className() |
| 462 | /// if classname is an empty string. |
| 463 | /// |
| 464 | /// If no object of the specified type is found, a NULL pointer is |
| 465 | /// returned. |
| 466 | |
| 467 | t = NULL__null; |
| 468 | |
| 469 | //map<const JObject*, string>::const_iterator iter = associated.begin(); |
| 470 | //for(; iter!=associated.end(); iter++){ |
| 471 | for( auto obj : associated ){ |
| 472 | t = dynamic_cast<const T*>(obj); |
| 473 | if(t!=NULL__null)return; |
| 474 | } |
| 475 | } |
| 476 | |
| 477 | //-------------------------- |
| 478 | // toStrings |
| 479 | //-------------------------- |
| 480 | inline void JObject::toStrings(vector<pair<string,string> > &items) const |
| 481 | { |
| 482 | /// Fill the given "items" vector with items representing the (important) |
| 483 | /// data members of this object. The structure of "items" is a vector |
| 484 | /// of pairs. The "first" element of the pair is the name of the item |
| 485 | /// as it should be displayed when dumping the item to the screen. For |
| 486 | /// example, one may wish to include units using a string like "r (cm)". |
| 487 | /// The "second" element of the pair is a formatted string containing the |
| 488 | /// value as it should be displayed. |
| 489 | /// |
| 490 | /// To facilitate this, the AddString() method exists which allows |
| 491 | /// items to be added with the desired formatting using a single line. |
| 492 | /// |
| 493 | /// This is a virtual method that is expected (but not required) |
| 494 | /// to be implemented by all classes that inherit from JObject. |
| 495 | |
| 496 | AddString(items, "JObject", "0x%08x", (unsigned long)this); |
| 497 | } |
| 498 | |
| 499 | //-------------------------- |
| 500 | // AddString |
| 501 | //-------------------------- |
| 502 | template<typename T> |
| 503 | void JObject::AddString(vector<pair<string,string> > &items, const char *name, const char *format, const T &val) const |
| 504 | { |
| 505 | /// Write the given value (val) to a string using the sprintf style formatting |
| 506 | /// string (format) and add it to the given vector (items) with the column |
| 507 | /// name "name". This is intended for use in the toStrings() method of |
| 508 | /// classes that inherit from JObject. |
| 509 | /// |
| 510 | /// The append_type flag provides a facility for recording the data type |
| 511 | /// and value with default formatting into items. This can be used |
| 512 | /// by a generic convertor (not part of JANA) to auto-generate a |
| 513 | /// representation of this object for use in some other persistence |
| 514 | /// package (e.g. ROOT files). |
| 515 | /// |
| 516 | /// If the append_types flag is set then the data type of "val" is |
| 517 | /// automatically appended with a colon (:) separator to the |
| 518 | /// name (first) part of the pair. In addition, "val" is converted |
| 519 | /// using stringstream and appended as well, also with a colon (:) |
| 520 | /// separator. For example, if the value of name passed in is "px" |
| 521 | /// and T is of type double, then the first member of the pair |
| 522 | /// appended to items will be something like "px:double:1.23784" |
| 523 | /// which can be decifered later to get the name, type, and value |
| 524 | /// of the data member. |
| 525 | /// |
| 526 | /// By default, the append_types flag is not set and the name part |
| 527 | /// of the pair is a straight copy of the name argument that is |
| 528 | /// passed in. |
| 529 | |
| 530 | char str[256]; |
| 531 | sprintf(str, format, val); |
| 532 | |
| 533 | stringstream ss; |
| 534 | ss<<name; |
| 535 | if(append_types){ |
| 536 | if(typeid(T)==typeid(int)){ |
| 537 | ss<<":int:"<<val; |
| 538 | }else if(typeid(T)==typeid(int32_t)){ |
| 539 | ss<<":int:"<<val; |
| 540 | }else if(typeid(T)==typeid(unsigned int)){ |
| 541 | ss<<":uint:"<<val; |
| 542 | }else if(typeid(T)==typeid(uint32_t)){ |
| 543 | ss<<":uint:"<<val; |
| 544 | }else if(typeid(T)==typeid(long)){ |
| 545 | ss<<":long:"<<val; |
| 546 | }else if(typeid(T)==typeid(int64_t)){ |
| 547 | ss<<":long:"<<val; |
| 548 | }else if(typeid(T)==typeid(unsigned long)){ |
| 549 | ss<<":ulong:"<<val; |
| 550 | }else if(typeid(T)==typeid(uint64_t)){ |
| 551 | ss<<":ulong:"<<val; |
| 552 | }else if(typeid(T)==typeid(short)){ |
| 553 | ss<<":short:"<<val; |
| 554 | }else if(typeid(T)==typeid(int16_t)){ |
| 555 | ss<<":short:"<<val; |
| 556 | }else if(typeid(T)==typeid(unsigned short)){ |
| 557 | ss<<":ushort:"<<val; |
| 558 | }else if(typeid(T)==typeid(uint16_t)){ |
| 559 | ss<<":ushort:"<<val; |
| 560 | }else if(typeid(T)==typeid(float)){ |
| 561 | ss<<":float:"<<val; |
| 562 | }else if(typeid(T)==typeid(double)){ |
| 563 | ss<<":double:"<<val; |
| 564 | }else if(typeid(T)==typeid(string)){ |
| 565 | ss<<":string:"<<val; |
| 566 | }else if(typeid(T)==typeid(const char*)){ |
| 567 | ss<<":string:"<<val; |
| 568 | }else if(typeid(T)==typeid(char*)){ |
| 569 | ss<<":string:"<<val; |
| 570 | }else{ |
| 571 | ss<<":unknown:"<<str; |
| 572 | } |
| 573 | } |
| 574 | |
| 575 | pair<string, string> item; |
| 576 | item.first = ss.str(); |
| 577 | item.second = string(str); |
| 578 | items.push_back(item); |
| 579 | } |
| 580 | |
| 581 | #endif // __CINT__ __CLING__ |
| 582 | |
| 583 | |
| 584 | } // Close JANA namespace |
| 585 | |
| 586 | #endif // _JObject_ |
| 587 |