Commit f0fa644a075ddc27bf11870bb0a426f6ebf1537d

Authored by Marcel Hohmann
1 parent 20941aa3

working clustering system

DualTrackTagger/include/DualTrackTaggerModule.h
... ... @@ -49,7 +49,7 @@ namespace Belle2 {
49 49 std::vector<Particle*> m_ClusterParticles;
50 50 };
51 51  
52   - /** This module selects and tags possible dual tracks as per the cuts outlined
  52 + /** This module selects and tags possible curl tracks as per the cuts outlined
53 53 * in Belle Note #1079 section 4.1.
54 54 * Currently works only for Belle data (with BelleFlag = True) converted via
55 55 * b2bii.
... ... @@ -83,7 +83,7 @@ namespace Belle2 {
83 83  
84 84 /** Loops through all particles in event adding those below a Pt threshold
85 85 to \code m_candidateParticles \endcode */
86   - virtual void identifyPotentialDuplicateTracks(Particle* p);
  86 + virtual void preSelectPotentialDuplicateTracks(Particle* p);
87 87  
88 88  
89 89 /** Loops through candidate tracks to compare and tag tracks if required */
... ... @@ -100,6 +100,8 @@ namespace Belle2 {
100 100  
101 101 virtual void addToClusters(Particle *p);
102 102  
  103 + virtual void tagClusterParticles();
  104 +
103 105 // ---------------- Functions related to Histos ------------
104 106  
105 107 /** define the histograms - only if MCFlag = true */
... ...
DualTrackTagger/src/DualTrackTaggerModule.cc
... ... @@ -64,20 +64,31 @@
64 64 if (!particlelist) {
65 65 B2ERROR("ParticleList " << iList << " not found");
66 66 continue;
67   - } else {
68   - if (particlelist->getListSize() == 0) {
69   - continue;
70   - }
71   - for (unsigned int i = 0; i < particlelist->getListSize(); ++i) {
72   - Particle* iParticle = particlelist->getParticle(i);
73   - m_allParticlesInList.push_back(iParticle);
74   - identifyPotentialDuplicateTracks(iParticle);
75   - } // loop through particles
76   - loopTracks();
77   - fill_h_genParticleIndexMulti();
78   - } // particle list exists
79   - } // loop list of particle lists
  67 + }
  68 + if (particlelist->getListSize() == 0) {
  69 + continue;
  70 + }
  71 + for (unsigned int i = 0; i < particlelist->getListSize(); ++i) {
  72 + Particle* iParticle = particlelist->getParticle(i);
  73 + m_allParticlesInList.push_back(iParticle);
  74 + preSelectPotentialDuplicateTracks(iParticle); // preselection
  75 + } // loop through particles
  76 +
  77 + for (Particle* p : m_candidateParticles) {
  78 + addToClusters(p);
  79 + }
  80 + tagClusterParticles();
  81 + //loopTracks();
  82 + fill_h_genParticleIndexMulti();
  83 +
  84 + for (DualTrackTaggerCluster *c : m_clusters) {
  85 + c -> saveClusterSizeToParticleExtraInfo();
  86 + }
  87 +
  88 + m_candidateParticles.clear();
  89 + m_clusters.clear();
80 90  
  91 + } // loop list of particle lists
81 92 } // DualTrackTaggerModule::event()
82 93  
83 94 void DualTrackTaggerModule::endRun()
... ... @@ -87,21 +98,25 @@
87 98 void DualTrackTaggerModule::terminate()
88 99 {
89 100 finalizeHistograms();
90   - for (DualTrackTaggerCluster *c : m_clusters) {
91   - c -> saveClusterSizeToParticleExtraInfo();
92   - }
93 101 //std::cout << "FINAL STATISTICS:" << total_particles << " | " << lowPt_particles << std::endl;
94 102 //std::cout << "FINAL STATISTICS:" << total_events << std::endl;
95 103 //std::cout << "FINAL STATISTICS:" << (float)total_particles / total_events << " | " << (float)lowPt_particles / total_events << std::endl;
96 104 }
97 105  
98   - void DualTrackTaggerModule::identifyPotentialDuplicateTracks(Particle* p)
  106 + void DualTrackTaggerModule::preSelectPotentialDuplicateTracks(Particle* p)
99 107 {
100 108 //Initiate all particles as not being dual tracks
101   - if (!p->hasExtraInfo("DualTrack")) {
102   - p -> addExtraInfo("DualTrack",0);
  109 + if (!p->hasExtraInfo("CurlTrack")) {
  110 + p -> addExtraInfo("CurlTrack",0);
103 111 } else {
104   - p -> setExtraInfo("DualTrack",0); // This should not be needed
  112 + B2ERROR("ERROR: particle already has CurlTrack extra info");
  113 + //p -> setExtraInfo("CurlTrack",0); // This should not be needed
  114 + }
  115 +
  116 + if (!p->hasExtraInfo("ClusterSize")) {
  117 + p -> addExtraInfo("ClusterSize",0);
  118 + } else {
  119 + B2ERROR("ERROR: particle already has ClusterSize extra info");
105 120 }
106 121  
107 122 // select cut for Belle or Belle 2
... ... @@ -186,6 +201,7 @@
186 201 double pTp1 = Variable::particlePt(p1);
187 202 double pTp2 = Variable::particlePt(p2);
188 203 double phi = p1->getMomentum().Angle(p2->getMomentum());
  204 + int absPDGp1 = abs(p1 -> getPDGcode());
189 205  
190 206 // get BDT response
191 207 double MVAresponseValue = Belle2ClassifierResponse -> GetMvaValue( {magDiffP, chargeMult, phi, pTp1, pTp2} );
... ... @@ -203,7 +219,7 @@
203 219 void DualTrackTaggerModule::tagDuplicate(Particle* p1, Particle* p2)
204 220 {
205 221 //prevent double tagging ---------- Allow double tagging for cases where there are multiple (>2) tracks. At least 1 will always be tagged as a normal track.
206   - //if (p1 -> getExtraInfo("DualTrack") == 1 || p2 -> getExtraInfo("DualTrack") == 1) {
  222 + //if (p1 -> getExtraInfo("CurlTrack") == 1 || p2 -> getExtraInfo("CurlTrack") == 1) {
207 223 // return;
208 224 //}
209 225  
... ... @@ -212,9 +228,9 @@
212 228 double p2chi2 = chi2(p2);
213 229  
214 230 if (p1chi2 > p2chi2) {
215   - p1 -> setExtraInfo("DualTrack",1);
  231 + p1 -> setExtraInfo("CurlTrack",1);
216 232 } else {
217   - p2 -> setExtraInfo("DualTrack",1);
  233 + p2 -> setExtraInfo("CurlTrack",1);
218 234 }
219 235 }
220 236  
... ... @@ -354,22 +370,47 @@
354 370 void DualTrackTaggerModule::addToClusters(Particle* p)
355 371 {
356 372 bool addedToACluster = false;
357   -
358 373 for (DualTrackTaggerCluster *cluster : m_clusters) {
359 374 if (matchesAllParticlesInCluster(p, cluster)) {
360 375 cluster->addToCluster(p);
  376 + //std::cout << "Added to existing cluster of size " << cluster -> size() << std::endl;
361 377 addedToACluster = true;
362 378 break;
363 379 }
364 380 }
365 381  
366 382 if (!addedToACluster) {
367   - DualTrackTaggerCluster *tempCluster;
  383 + DualTrackTaggerCluster *tempCluster = new DualTrackTaggerCluster();
368 384 tempCluster->addToCluster(p);
  385 + //std::cout << tempCluster -> size() << std::endl;
369 386 m_clusters.push_back(tempCluster);
370 387 }
371 388 }
372 389  
  390 +void DualTrackTaggerModule::tagClusterParticles()
  391 +{
  392 + for (DualTrackTaggerCluster *c : m_clusters) {
  393 + // find best chi2
  394 + int posLowestChi2 = 0;
  395 + double LowestChi2 = 0;
  396 +
  397 + for (int i = 0; i < c -> size(); i++) {
  398 + if (chi2(c -> getParticle(i)) < LowestChi2) {
  399 + LowestChi2 = chi2(c -> getParticle(i));
  400 + posLowestChi2 = i;
  401 + }
  402 + } // particles within cluster
  403 +
  404 + for (int i = 0; i < c -> size(); i++) {
  405 + if (i == posLowestChi2) {
  406 + c -> getParticle(i) -> setExtraInfo("CurlTrack",0);
  407 + } else {
  408 + c -> getParticle(i) -> setExtraInfo("CurlTrack",1);
  409 + }
  410 + } // particles within cluster
  411 + } // clusters
  412 +}
  413 +
373 414  
374 415 bool DualTrackTaggerModule::matchesAllParticlesInCluster(Particle* p, DualTrackTaggerCluster *c)
375 416 {
... ... @@ -405,10 +446,11 @@ void DualTrackTaggerCluster::saveClusterSizeToParticleExtraInfo()
405 446 {
406 447 int clusterSize = m_ClusterParticles.size();
407 448 for (Particle *p : m_ClusterParticles) {
408   - if (!p->hasExtraInfo("ClusterSize")) {
409   - p -> addExtraInfo("ClusterSize",clusterSize);
410   - } else { // in case the particle is in two clusters
411   - p -> setExtraInfo("ClusterSize", 0.5*(clusterSize + p->getExtraInfo("ClusterSize")));
412   - }
  449 + p -> setExtraInfo("ClusterSize",clusterSize);
  450 + //if (!p->hasExtraInfo("ClusterSize")) {
  451 + //} else { // in case the particle is in two clusters (shouldnt be possible)
  452 + //std:: cout << "ERROR: particle registered in two clusters" << std::endl;
  453 + //p -> setExtraInfo("ClusterSize", 0.5*(clusterSize + p->getExtraInfo("ClusterSize")));
  454 + //}
413 455 }
414 456 }
... ...