ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/HbbAnalysis/src/JetFlavour.cc
(Generate patch)

Comparing UserCode/HbbAnalysis/src/JetFlavour.cc (file contents):
Revision 1.1 by amagnan, Thu May 28 15:12:45 2009 UTC vs.
Revision 1.9 by amagnan, Thu Mar 17 12:47:39 2011 UTC

# Line 49 | Line 49 | namespace HbbAnalysis {
49      flav_.clear();
50      eid_.clear();
51      muid_.clear();
52 <    partons_.reserve(10);
53 <    eid_.reserve(10);
54 <    muid_.reserve(10);
55 <    flav_.reserve(10);
52 >    partons_.reserve(2);
53 >    eid_.reserve(2);
54 >    muid_.reserve(2);
55 >    flav_.reserve(2);
56  
57      unsigned int n_quarks = 0;
58      unsigned int n_antiquarks = 0;
# Line 64 | Line 64 | namespace HbbAnalysis {
64  
65      if (debug_) std::cout << "* nMC = " << lNMC << std::endl;
66  
67 <    //loop on MC part to find b and bbar initial particles.
67 >    //loop on MC part to find any initial partons.
68      for (unsigned int imc(0); imc < lNMC; imc++){//loop on MC particles
69        const reco::GenParticle & lPart = (*aGenParticles).at(imc);
70        //break loop if W events
# Line 76 | Line 76 | namespace HbbAnalysis {
76               abs(lPart.pdgId()) == 16 ) ) {
77          isW = true;
78          nWEvts_++;
79 <        if (!isZ) break;
79 >        //if (!isZ) break;
80        }
81        //record when Z event
82        if ( lPart.status() == 3 && lPart.numberOfMothers() > 1 &&
# Line 90 | Line 90 | namespace HbbAnalysis {
90        //looking for the final-states partons, which will have 2 parents being the 2 initial states partons.
91        if (lPart.status() == 3 && lPart.numberOfMothers() > 1 &&
92            ( ( flavour_ == 45 && (abs(lPart.pdgId()) == 4 || abs(lPart.pdgId()) == 5) ) ||
93 <            ( flavour_ > 3 && flavour_ < 22 && abs(lPart.pdgId()) == flavour_) ||
94 <            ( flavour_ == 3 && ( (abs(lPart.pdgId()) <= flavour_ ||
93 >            ( flavour_ > 3 && flavour_ < 22 && static_cast<unsigned int>(abs(lPart.pdgId())) == flavour_) ||
94 >            ( flavour_ == 3 && ( (static_cast<unsigned int>(abs(lPart.pdgId())) <= flavour_ ||
95                                    lPart.pdgId() == 21)
96                                   )
97                )
98              )
99            ) {
100 <        if (debug_ > 0) {
100 >        if (debug_ > 1) {
101            std::cout << "** initial parton #" << imc ;
102            printParticle(lPart);
103          }
# Line 106 | Line 106 | namespace HbbAnalysis {
106          else n_antiquarks++;
107          p_partonFlavour->Fill(lPart.pdgId());
108        }
109 <      else if ((debug_ > 0 && imc < 15) || debug_ > 1) {
109 >      else if ((debug_ > 1 && imc < 15) || debug_ > 2) {
110          std::cout << "** not selected: #" << imc;
111          printParticle(lPart);
112        }
# Line 114 | Line 114 | namespace HbbAnalysis {
114      }//loop on MC particles
115  
116      //if (isW == isZ) std::cout << "*** WARNING *** isW = " << isW << ", isZ = " << isZ << std::endl;
117 <    if ((flavour_ > 3 && flavour_ != 21) && !isZ) return 0;
117 >    //if ((flavour_ > 3 && flavour_ != 21) && !isZ) return 0;
118  
119      if (debug_) {
120        std::cout << "nQuarks = " << n_quarks << ", n_antiQuarks = " << n_antiquarks << std::endl;
# Line 136 | Line 136 | namespace HbbAnalysis {
136      nTotal_ += nPartons;
137      assert (nPartons == n_quarks+n_antiquarks);
138  
139 <    if (flavour_ > 3) {
140 <      assert(n_quarks >= 1 && n_antiquarks >= 1);
141 <      assert (nPartons == 2);
142 <    }
139 > //    if (flavour_ > 3) {
140 > //       assert(n_quarks >= 1 && n_antiquarks >= 1);
141 > //       assert (nPartons == 2);
142 > //     }
143  
144      p_nQuarks->Fill(n_quarks);
145      p_nAntiQuarks->Fill(n_antiquarks);
# Line 158 | Line 158 | namespace HbbAnalysis {
158      unsigned int *eid  = new unsigned int[nPartons];
159      unsigned int *muid = new unsigned int[nPartons];
160  
161 –
161      for (unsigned int ip(0); ip<nPartons; ip++){//loop on partons
162        nDaugh[ip] = 0;
163        flav[ip] = 1;
# Line 169 | Line 168 | namespace HbbAnalysis {
168        eid[ip] = 0;
169        muid[ip] = 0;
170        daughters[ip].clear();
171 <      daughters[ip].reserve(250);
171 >      daughters[ip].reserve(500);
172        parents[ip].clear();
173 <      parents[ip].reserve(250);
173 >      parents[ip].reserve(500);
174        parents[ip].push_back(partons_.at(ip));
175  
176        for (unsigned int pmc = 0; pmc < parents[ip].size(); pmc++){//loop "interactively" on parents
177          unsigned int lParIndex = parents[ip].at(pmc);
178 <        if (debug_ > 2) std::cout << "** Looking for parent : #" << lParIndex << " (id " <<  (*aGenParticles).at(lParIndex).pdgId() << "), vector size = " <<  parents[ip].size() << std::endl;
178 >        const reco::GenParticle & ppar = (*aGenParticles).at(lParIndex);
179 >        if (debug_ > 2) std::cout << "** Looking for parent : #" << lParIndex << " (id " <<  ppar.pdgId() << "), vector size = " <<  parents[ip].size() << std::endl;
180  
181          //counter to check a match is found
182          unsigned int lNNotFound = 0;
# Line 187 | Line 187 | namespace HbbAnalysis {
187  
188            for (unsigned int im(0); im<pgen.numberOfMothers(); im++){//loop on parents
189              //recursive loop to find all parents
190 <            if ( compareParticles(pgen.mother(im),(*aGenParticles).at(lParIndex)) ){//if par=parton
190 >            if ( compareParticles(pgen.mother(im),ppar)){//if par=parton
191                lNTot++;
192                if (pgen.status() != 1) {
193                  nDaugh[ip]++;
# Line 227 | Line 227 | namespace HbbAnalysis {
227          }//loop on parts
228          if (lNNotFound == lNTot ) {
229            std::cout << " --- WARNING: no particle found for parton #" << lParIndex ;
230 <          printParticle((*aGenParticles).at(lParIndex));
230 >          printParticle(ppar);
231          }
232  
233        }//loop on parents
# Line 264 | Line 264 | namespace HbbAnalysis {
264        }//for b and c jets
265  
266  
267 <      if (parents[ip].size() > 250) std::cout << "WARNING ! Initial default size of parents not sufficient = " << parents[ip].size() << std::endl;
268 <      if (daughters[ip].size() > 250) std::cout << "WARNING ! Initial default size of daughters not sufficient = " << daughters[ip].size() << std::endl;
267 >      if (parents[ip].size() > 500) std::cout << "WARNING ! Initial default size of parents not sufficient = " << parents[ip].size() << std::endl;
268 >      if (daughters[ip].size() > 500) std::cout << "WARNING ! Initial default size of daughters not sufficient = " << daughters[ip].size() << std::endl;
269  
270  
271      }//loop on partons
# Line 303 | Line 303 | namespace HbbAnalysis {
303  
304      }//for light jets
305  
306 +
307      delete[] nDaugh;
308      delete[] parents;
309      delete[] daughters;
# Line 316 | Line 317 | namespace HbbAnalysis {
317    }//fillPartons
318  
319  
320 <  std::pair<int,unsigned int> JetFlavour::partonMatchingGenJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, double aDR){//parton match with gen jet
320 >  const std::pair<int,unsigned int> JetFlavour::partonMatchingGenJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, const double aDR) const{//parton match with gen jet
321  
322      int aIndex = -1;
323      double deltaRMin = 1000;
# Line 336 | Line 337 | namespace HbbAnalysis {
337  
338    }//parton match with gen jet
339  
340 <  unsigned int JetFlavour::leptonMatchingGenJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, double aDR){//lepton match with gen jet
340 >  const unsigned int JetFlavour::leptonMatchingGenJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, const double aDR) const{//lepton match with gen jet
341  
342      int aIndex = -1;
343      double deltaRMin = 1000;
# Line 372 | Line 373 | namespace HbbAnalysis {
373  
374    }//lepton match with gen jet
375  
376 <  unsigned int JetFlavour::partonMatchingRecoJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, double aDR){//parton match with reco jet
376 >  const unsigned int JetFlavour::partonMatchingRecoJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, const double aDR) const{//parton match with reco jet
377  
378      int aIndex = -1;
379      double deltaRMin = 1000;
# Line 390 | Line 391 | namespace HbbAnalysis {
391  
392    }//parton match with reco jet
393  
394 <  unsigned int JetFlavour::leptonMatchingRecoJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, double aDR){//lepton match with reco jet
394 >  const unsigned int JetFlavour::leptonMatchingRecoJet(const pat::Jet & aJet, const edm::Handle<reco::GenParticleCollection> & aGenParticles, const double aDR) const{//lepton match with reco jet
395  
396      int aIndex = -1;
397      double deltaRMin = 1000;
# Line 454 | Line 455 | namespace HbbAnalysis {
455      
456    }
457  
458 <  void JetFlavour::printSummary() {
458 >  void JetFlavour::printSummary() const{
459  
460      std::cout <<
461        "=====================================================================" << std::endl <<

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines