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

Comparing UserCode/algomez/GenAnalyzer/src/GenAnalyzer.cc (file contents):
Revision 1.1 by algomez, Mon Jun 3 20:36:36 2013 UTC vs.
Revision 1.2 by algomez, Wed Jun 19 02:42:01 2013 UTC

# Line 67 | Line 67 | class GenAnalyzer : public edm::EDAnalyz
67          edm::InputTag src_;
68          double stop1Mass_;
69          double stop2Mass_;
70 +        double st1decay_;
71  
72          TH1D * histo;
73          // Higgs
# Line 148 | Line 149 | class GenAnalyzer : public edm::EDAnalyz
149   GenAnalyzer::GenAnalyzer(const edm::ParameterSet& iConfig) :
150          src_( iConfig.getParameter<edm::InputTag>( "src" ) ),                   // Obtain inputs
151          stop1Mass_( iConfig.getParameter<double>( "stop1Mass" ) ),
152 <        stop2Mass_( iConfig.getParameter<double>( "stop2Mass" ) )
152 >        stop2Mass_( iConfig.getParameter<double>( "stop2Mass" ) ),
153 >        st1decay_( iConfig.getParameter<double>( "st1decay" ) )
154   {
155          //now do what ever initialization is needed
156          edm::Service<TFileService> fs;                                          // Output File
# Line 264 | Line 266 | void GenAnalyzer::analyze(const edm::Eve
266  
267          // Begin Loop for GenParticles
268          for(unsigned int i = 0; i < particles->size(); ++ i) {
269 +
270 +                if( p[i].status()!= 3 ) continue;                               // Make sure only "final" particles are present
271                  
272                  const Candidate * mom = p[i].mother();                          // call mother particle
273 <                int idAbs = abs( p[i].pdgId() );                                // call pdg Id
270 <                if( idAbs == 5 ) histo->Fill( p[i].pt() );                      // test
273 >                if( p[i].pdgId() == 5 ) histo->Fill( p[i].pt() );                       // test
274  
275  
276                  //////////// True Mass plots
277 <                if( idAbs == 25 ){
277 >                if( p[i].pdgId() == 25 ){
278                          higgs_num++;                                            // Check number of Higgs
279                          h_higgsTrue_mass->Fill(p[i].mass());                    // Higgs true mass
280                  }
281 <                if( idAbs == 1000002 ) h_stop1True_mass->Fill(p[i].mass());     // Stop1 true mass
282 <                if( idAbs == 1000004 ) h_stop2True_mass->Fill(p[i].mass());     // Stop2 true mass
281 >                if( p[i].pdgId() == 1000002 ) h_stop1True_mass->Fill(p[i].mass());      // Stop1 true mass
282 >                if( p[i].pdgId() == 1000004 ) h_stop2True_mass->Fill(p[i].mass());      // Stop2 true mass
283                  ////////////
284                  
282 –                if( p[i].status()!= 3 ) continue;                               // Make sure only "final" particles are present
285  
286                  //////////// Storing TLorentz Vectors
287                  TLorentzVector tmpjet, tmpBjet, tmpjetsB, tmpHiggs, tmpStop1B, tmpStop1jet;     // tmp TLorentzVectors
288                  if( mom ){
289                          // Jets without b
290 <                        if( ( idAbs == 1 ) | ( idAbs == 2 ) | ( idAbs == 3 ) | ( idAbs == 4) ){
290 >                        if( ( p[i].pdgId() == 1 ) || ( p[i].pdgId() == 2 ) || ( p[i].pdgId() == 3 ) || ( p[i].pdgId() == 4) ){
291                                  tmpjet.SetPxPyPzE(p[i].px(),p[i].py(),p[i].pz(),p[i].energy());
292                                  p4jet.push_back( tmpjet );
293                          }
294                          std::sort(p4jet.begin(), p4jet.end(), ComparePt);
295                          // Bjets
296 <                        if( idAbs == 5 ){
296 >                        if( p[i].pdgId() == 5 ){
297                                  tmpBjet.SetPxPyPzE(p[i].px(),p[i].py(),p[i].pz(),p[i].energy());
298                                  p4Bjet.push_back( tmpBjet );
299                          }
300                          std::sort(p4Bjet.begin(), p4Bjet.end(), ComparePt);
301                          // Jets with b
302 <                        if( ( idAbs == 1 ) | ( idAbs == 2 ) | ( idAbs == 3 ) | ( idAbs == 4) | ( idAbs == 5 )){
302 >                        if( ( p[i].pdgId() == 1 ) || ( p[i].pdgId() == 2 ) || ( p[i].pdgId() == 3 ) || ( p[i].pdgId() == 4) || ( p[i].pdgId() == 5 )){
303                                  tmpBjet.SetPxPyPzE(p[i].px(),p[i].py(),p[i].pz(),p[i].energy());
304                                  p4jetsB.push_back( tmpBjet );
305                          }
306                          std::sort(p4jetsB.begin(), p4jetsB.end(), ComparePt);
307  
308                          // Higgs
309 <                        if( ( idAbs == 5 ) & ( mom->pdgId() == 25 ) ){
309 >                        if( ( p[i].pdgId() == 5 ) && ( mom->pdgId() == 25 ) ){
310                                  tmpHiggs.SetPxPyPzE(p[i].px(),p[i].py(),p[i].pz(),p[i].energy());
311                                  p4HiggsB.push_back( tmpHiggs );
312                          }
313                          std::sort(p4HiggsB.begin(), p4HiggsB.end(), ComparePt);
314                          // Stop1
315                          if( abs( mom->pdgId() )== 1000002 ){
316 <                                if( ( idAbs == 1 ) | ( idAbs == 2 ) | ( idAbs == 3 ) | ( idAbs == 4) ){
315 <                                        tmpStop1B.SetPxPyPzE(p[i].px(),p[i].py(),p[i].pz(),p[i].energy());
316 <                                        p4Stop1B.push_back( tmpStop1B );
317 <                                        p4Stop1jetsB.push_back( tmpStop1B );
318 <                                }
319 <                                else {
316 >                                if( ( p[i].pdgId() == 1 ) || ( p[i].pdgId() == 2 ) || ( p[i].pdgId() == 3 ) || ( p[i].pdgId() == 4) ){
317                                          tmpStop1jet.SetPxPyPzE(p[i].px(),p[i].py(),p[i].pz(),p[i].energy());
318                                          p4Stop1jet.push_back( tmpStop1jet );
319                                          p4Stop1jetsB.push_back( tmpStop1jet );
320                                  }
321 +                                else {
322 +                                        tmpStop1B.SetPxPyPzE(p[i].px(),p[i].py(),p[i].pz(),p[i].energy());
323 +                                        p4Stop1B.push_back( tmpStop1B );
324 +                                        p4Stop1jetsB.push_back( tmpStop1B );
325 +                                }
326                          }
327                          std::sort(p4Stop1B.begin(), p4Stop1B.end(), ComparePt);
328                          std::sort(p4Stop1jet.begin(), p4Stop1jet.end(), ComparePt);
# Line 328 | Line 330 | void GenAnalyzer::analyze(const edm::Eve
330                  }
331          }  // end GenParticles loop
332  
333 +
334          // Test number of true Higgs
335          h_higgs_num->Fill( higgs_num );
336  
337 +        //cout << "run:lumi:event = " << iEvent.id().run() << ":" << iEvent.id().luminosityBlock() << ":" << iEvent.id().event() <<endl;
338  
339          ///// Higgs reconstructed mass
340          for(unsigned int j = 0; j < p4HiggsB.size()-1; ++j) {
# Line 349 | Line 353 | void GenAnalyzer::analyze(const edm::Eve
353                          }
354                  }
355          }
356 <        //////////////////////////////
356 >        ///////////////////////////// /
357  
358  
359          /////// Stop1 reconstructed mass
360 <        for(unsigned int j = 0; j < p4Stop1B.size(); ++j) {
361 <                for(unsigned int k= 0; k < p4Stop1jet.size(); ++k) {
362 <                        TLorentzVector p4CandidateStop1 = p4Stop1B[j] + p4Stop1jet[k];          // Higgs TLorentzVector
363 <                        Double_t massStop1 = stop1Mass_ ;                                               // nominal mass value
364 <                        Double_t diffmassStop1 = p4CandidateStop1.M() - massStop1;              // mass difference
365 <                        if( abs(diffmassStop1) < 1.0 ){                                         // diffmass presicion
366 <                                h_stop1_mass->Fill( p4CandidateStop1.M() );                     // Stop1 mass
367 <                                h_stop1_pt->Fill( p4CandidateStop1.Pt() );                      // Stop1 pt
368 <                                p4Stop1.push_back( p4CandidateStop1 );
360 >        if( st1decay_==1){
361 >                for(unsigned int j = 0; j < p4Stop1B.size(); ++j) {
362 >                        for(unsigned int k= 0; k < p4Stop1jet.size(); ++k) {
363 >                                TLorentzVector p4CandidateStop1 = p4Stop1B[j] + p4Stop1jet[k];          // higgs tlorentzvector
364 >                                Double_t massStop1 = stop1Mass_ ;                                               // nominal mass value
365 >                                Double_t diffmassStop1 = p4CandidateStop1.M() - massStop1;              // mass difference
366 >                                if( abs(diffmassStop1) < 1.0 ){                                         // diffmass presicion
367 >                                        h_stop1_mass->Fill( p4CandidateStop1.M() );                     // stop1 mass
368 >                                        h_stop1_pt->Fill( p4CandidateStop1.Pt() );                      // stop1 pt
369 >                                        p4Stop1.push_back( p4CandidateStop1 );
370 >                                }
371 >                        }
372 >                }
373 >        } else {
374 >                for(unsigned int j = 0; j < p4Stop1jet.size(); ++j) {
375 >                        for(unsigned int k= 0; k < p4Stop1jet.size(); ++k) {
376 >                                if( j==k ) continue;
377 >                                TLorentzVector p4CandidateStop1 = p4Stop1jet[j] + p4Stop1jet[k];                // higgs tlorentzvector
378 >                                Double_t massStop1 = stop1Mass_ ;                                               // nominal mass value
379 >                                Double_t diffmassStop1 = p4CandidateStop1.M() - massStop1;              // mass difference
380 >                                if( abs(diffmassStop1) < 1.0 ){                                         // diffmass presicion
381 >                                        h_stop1_mass->Fill( p4CandidateStop1.M() );                     // stop1 mass
382 >                                        h_stop1_pt->Fill( p4CandidateStop1.Pt() );                      // stop1 pt
383 >                                        p4Stop1.push_back( p4CandidateStop1 );
384 >                                }
385                          }
386                  }
387          }
388 <        //////////////////////////////////
388 >        ///////////////////////////////// /
389          
390          // Stop2 reconstructed mass
391          for(unsigned int j = 0; j < p4Higgs.size(); ++j) {
# Line 377 | Line 397 | void GenAnalyzer::analyze(const edm::Eve
397                                  h_stop2_mass->Fill( p4CandidateStop2.M() );                     // Stop1 mass
398                          }
399                  }
400 <        }
400 >        }//
401  
402          // Jets + B
403          Double_t Ht = 0;
# Line 405 | Line 425 | void GenAnalyzer::analyze(const edm::Eve
425                  for(unsigned int k= 0; k < p4Stop1jet.size(); ++k) {
426                          if( p4Stop1jet[k].Pt() == p4jetsB[j].Pt() ) h_jetfromStop1->Fill( j+1 );
427                                  //cout << k << " "<< j << " " <<p4Stop1jet[k].Pt() << " " << p4jetsB[j].Pt() << endl;
428 <                }
428 >                }/*
429                  if( p4Stop1B[0].Pt() == p4jetsB[j].Pt() ) h_Bjet1fromStop1->Fill( j+1 );
430 <                if( p4Stop1B[1].Pt() == p4jetsB[j].Pt() ) h_Bjet2fromStop1->Fill( j+1 );
430 >                if( p4Stop1B[1].Pt() == p4jetsB[j].Pt() ) h_Bjet2fromStop1->Fill( j+1 );*/
431                  for(unsigned int kk= 0; kk < p4Stop1B.size(); ++kk) {
432                          if( p4Stop1B[kk].Pt() == p4jetsB[j].Pt() ) h_BjetfromStop1->Fill( j+1 );
433                                  //cout << kk << " "<< j << " " <<p4Stop1B[kk].Pt() << " " << p4jetsB[j].Pt() << endl;
# Line 422 | Line 442 | void GenAnalyzer::analyze(const edm::Eve
442                  }
443          }
444          h_jetsB_ht->Fill( Ht );
445 <        /////////////////////////////
445 >        ////////////////////////////
446  
447          // Bjets
448          for(unsigned int j = 0; j < p4Bjet.size(); ++j) {
# Line 433 | Line 453 | void GenAnalyzer::analyze(const edm::Eve
453                  h_Bjet1Bjet2_cosDeltaPhi->Fill( cos( p4Bjet[0].DeltaPhi( p4Bjet[1] ) ) );
454                  h_Bjet1Bjet2_deltaR->Fill( p4Bjet[0].DeltaR( p4Bjet[1] ) );
455          }
456 <        /////////////////////////////
456 >        //////////////////////////// /
457  
458          // Jets
459          for(unsigned int j = 0; j < p4jet.size(); ++j) {
460                  h_jet_pt->Fill( p4jet[j].Pt() );
461                  h_jet1_pt->Fill( p4jet[0].Pt() );
462                  h_jet2_pt->Fill( p4jet[1].Pt() );
443 –                h_jet_num->Fill( p4jet.size() );
463          }
464 +                h_jet_num->Fill( p4jet.size() );
465          /////////////////////////////
466  
467   }

Diff Legend

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