Reference Manual

Go to the documentation of this file.
00001 
00002 #include "Validation/EventGenerator/interface/VVVValidation.h"
00003 #include "Validation/EventGenerator/interface/HepMCValidationHelper.h"
00004 
00005 
00006 #include "CLHEP/Units/defs.h"
00007 #include "CLHEP/Units/PhysicalConstants.h"
00008 
00009 
00010 using namespace edm;
00011 using namespace std;
00012 
00013 VVVValidation::VVVValidation(const edm::ParameterSet& iPSet): 
00014   _wmanager(iPSet),
00015   hepmcCollection_(iPSet.getParameter<edm::InputTag>("hepmcCollection")),
00016   genparticleCollection_(iPSet.getParameter<edm::InputTag>("genparticleCollection")),
00017   genjetCollection_(iPSet.getParameter<edm::InputTag>("genjetsCollection")),
00018   matchPr_(iPSet.getParameter<double>("matchingPrecision")),
00019   _lepStatus(iPSet.getParameter<double>("lepStatus")),
00020   verbosity_(iPSet.getUntrackedParameter<unsigned int>("verbosity",0))
00021 {    
00022   dbe = 0;
00023   dbe = edm::Service<DQMStore>().operator->();
00024 }
00025 
00026 VVVValidation::~VVVValidation() {}
00027 
00028 void VVVValidation::beginJob()
00029 {
00030 
00031 
00032 
00033   if(dbe){
00034 
00035     dbe->setCurrentFolder("Generator/VVVJets");
00036         
00037     
00038     TH1::SetDefaultSumw2();
00039     nEvt = dbe->book1D("nEvt", "n analyzed Events", 1, 0., 1.);
00040 
00041     genJetMult = dbe->book1D("genJetMult", "GenJet multiplicity", 50, 0, 50);
00042     genJetEnergy = dbe->book1D("genJetEnergy", "Log10(GenJet energy)", 60, -1, 5);
00043     genJetPt = dbe->book1D("genJetPt", "Log10(GenJet pt)", 60, -1, 5);
00044     genJetEta = dbe->book1D("genJetEta", "GenJet eta", 220, -11, 11);
00045     genJetPhi = dbe->book1D("genJetPhi", "GenJet phi", 360, -180, 180);
00046     genJetDeltaEtaMin = dbe->book1D("genJetDeltaEtaMin", "GenJet minimum rapidity gap", 30, 0, 30);
00047     
00048     genJetPto1 = dbe->book1D("genJetPto1", "GenJet multiplicity above 1 GeV", 50, 0, 50);
00049     genJetPto30 = dbe->book1D("genJetPto30", "GenJet multiplicity above 30 GeV", 10, -0.5, 9.5);
00050     genJetPto50 = dbe->book1D("genJetPto50", "GenJet multiplicity above 50 GeV", 10, -0.5, 9.5);
00051     genJetPto100 = dbe->book1D("genJetPto100", "GenJet multiplicity above 100 GeV", 10, -0.5, 9.5);
00052     genJetCentral = dbe->book1D("genJetCentral", "GenJet multiplicity |eta|.lt.2.5", 10, -0.5, 9.5);
00053 
00054     WW_TwoJEt_JetsM= dbe->book1D("WW_TwoJEt_JetsM", "2 jet invariant mass in WW 2 J events", 200, 0, 2000);
00055     h_l_jet_eta = dbe->book1D("h_l_jet_eta", "leading jet eta",  30, -5.2, 5.2);
00056     h_l_jet_pt = dbe->book1D("h_l_jet_pt", "leading jet pt", 50, 0, 300);
00057     h_sl_jet_eta = dbe->book1D("h_sl_jet_eta", "leading jet eta",  30, -5.2, 5.2);
00058     h_sl_jet_pt = dbe->book1D("h_sl_jet_pt", "leading jet pt", 50, 0, 300);
00059     h_ssl_jet_eta = dbe->book1D("h_ssl_jet_eta", "leading jet eta",  30, -5.2, 5.2); 
00060     h_ssl_jet_pt = dbe->book1D("h_ssl_jet_pt", "leading jet pt", 50, 0, 300);
00061 
00062 
00063     genJetTotPt = dbe->book1D("genJetTotPt", "Log10(GenJet total pt)", 100, 0, 1000);
00064  
00065     h_dr = dbe->book1D("h_dr", "DeltaR between leptons and jets", 30, 0, 2);
00066     h_mWplus = dbe->book1D("h_mWplus", "M W^{+}", 50, 0, 200);
00067     h_phiWplus= dbe->book1D("h_phiWplus", "#phi W^{+}", 30, -3.5, 3.5);
00068     h_ptWplus = dbe->book1D("h_ptWplus", "P_{T} W^{+}", 50, 0, 200);
00069     h_yWplus = dbe->book1D("h_yWplus", "Y W^{+}", 30, -3.5, 3.5);
00070 
00071     h_mWminus= dbe->book1D("h_mWminus", "M W^{-}", 50, 0, 200);
00072     h_phiWminus= dbe->book1D("h_phiWminus", "#phi W^{-}", 30, -3.5, 3.5);
00073     h_ptWminus= dbe->book1D("h_ptWminus", "P_{T} W^{-}", 50, 0, 200);
00074     h_yWminus= dbe->book1D("h_yWminus", "Y W^{-}", 30, -3.5, 3.5);
00075 
00076     h_mZ= dbe->book1D("h_mZ", "M Z", 50, 0, 200);
00077     h_phiZ= dbe->book1D("h_phiZ", "#phi Z", 30, -3.5, 3.5);
00078     h_ptZ= dbe->book1D("h_ptZ", "P_{T} Z", 50, 0, 200);
00079     h_yZ= dbe->book1D("h_yZ", "Y Z", 30, -3.5, 3.5);
00080 
00081     h_mWplus_3b = dbe->book1D("h_mWplus_3b", "M W^{+}", 50, 0, 200);
00082     h_phiWplus_3b= dbe->book1D("h_phiWplus_3b", "#phi W^{+}", 30, -3.5, 3.5);
00083     h_ptWplus_3b = dbe->book1D("h_ptWplus_3b", "P_{T} W^{+}", 50, 0, 200);
00084     h_yWplus_3b = dbe->book1D("h_yWplus_3b", "Y W^{+}", 30, -3.5, 3.5);
00085 
00086     h_mWminus_3b= dbe->book1D("h_mWminus_3b", "M W^{-}", 50, 0, 200);
00087     h_phiWminus_3b= dbe->book1D("h_phiWminus_3b", "#phi W^{-}", 30, -3.5, 3.5);
00088     h_ptWminus_3b= dbe->book1D("h_ptWminus_3b", "P_{T} W^{-}", 50, 0, 200);
00089     h_yWminus_3b= dbe->book1D("h_yWminus_3b", "Y W^{-}", 30, -3.5, 3.5);
00090 
00091     h_mZ_3b= dbe->book1D("h_mZ_3b", "M Z", 50, 0, 200);
00092     h_phiZ_3b= dbe->book1D("h_phiZ_3b", "#phi Z", 30, -3.5, 3.5);
00093     h_ptZ_3b= dbe->book1D("h_ptZ_3b", "P_{T} Z", 50, 0, 200);
00094     h_yZ_3b= dbe->book1D("h_yZ_3b", "Y Z", 30, -3.5, 3.5);
00095 
00096     h_mWW= dbe->book1D("h_mWW", "M W^{-}W^{+}", 200, 0, 2000);
00097     h_phiWW= dbe->book1D("h_phiWW", "#phi W^{-}W^{+}", 30, -3.5, 3.5);
00098     h_ptWW= dbe->book1D("h_ptWW", "P_{T} W^{-}W^{+}", 50, 0, 2000);
00099     h_yWW =dbe->book1D("h_yWW", "Y W^{-}W^{+}", 30, -3.5, 3.5);
00100 
00101     h_mWZ= dbe->book1D("h_mWZ", "M W Z", 200, 0, 2000);
00102     h_phiWZ= dbe->book1D("h_phiWZ", "#phi W Z", 30, -3.5, 3.5);
00103     h_ptWZ= dbe->book1D("h_ptWZ", "P_{T} W Z", 50, 0, 2000);
00104     h_yWZ =dbe->book1D("h_yWZ", "Y W Z", 30, -3.5, 3.5);
00105 
00106     h_mZZ= dbe->book1D("h_mZZ", "M Z Z", 200, 0, 2000);
00107     h_phiZZ= dbe->book1D("h_phiZZ", "#phi Z Z", 30, -3.5, 3.5);
00108     h_ptZZ= dbe->book1D("h_ptZZ", "P_{T} Z Z", 50, 0, 2000);
00109     h_yZZ =dbe->book1D("h_yZZ", "Y Z Z", 30, -3.5, 3.5);
00110 
00111     h_mWWW= dbe->book1D("h_mWWW", "M W W W ", 200, 0, 2000);
00112     h_phiWWW= dbe->book1D("h_phiWWW", "#phi W W W", 30, -3.5, 3.5);
00113     h_ptWWW= dbe->book1D("h_ptWWW", "P_{T} W W W", 50, 0, 2000);
00114     h_yWWW =dbe->book1D("h_yWWW", "Y W W W", 30, -3.5, 3.5);
00115 
00116     h_mWWZ= dbe->book1D("h_mWWZ", "M W W Z ", 200, 0, 2000);
00117     h_phiWWZ= dbe->book1D("h_phiWWZ", "#phi W W Z", 30, -3.5, 3.5);
00118     h_ptWWZ= dbe->book1D("h_ptWWZ", "P_{T} W W Z", 50, 0, 2000);
00119     h_yWWZ =dbe->book1D("h_yWWZ", "Y W W Z", 30, -3.5, 3.5);
00120 
00121     h_mWZZ= dbe->book1D("h_mWZZ", "M W Z Z ", 200, 0, 2000);
00122     h_phiWZZ= dbe->book1D("h_phiWZZ", "#phi W Z Z", 30, -3.5, 3.5);
00123     h_ptWZZ= dbe->book1D("h_ptWZZ", "P_{T} W Z Z", 50, 0, 2000);
00124     h_yWZZ =dbe->book1D("h_yWZZ", "Y W Z Z", 30, -3.5, 3.5);
00125 
00126     h_mZZZ= dbe->book1D("h_mZZZ", "M Z Z Z ", 200, 0, 2000);
00127     h_phiZZZ= dbe->book1D("h_phiZZZ", "#phi Z Z Z", 30, -3.5, 3.5);
00128     h_ptZZZ= dbe->book1D("h_ptZZZ", "P_{T} Z Z Z", 50, 0, 2000);
00129     h_yZZZ =dbe->book1D("h_yZZZ", "Y Z Z Z", 30, -3.5, 3.5);
00130 
00131  
00132           
00133     leading_l_pt = dbe->book1D("leading_l_pt", "leading lepton pt", 50, 0, 200);
00134     subleading_l_pt = dbe->book1D("subleading_l_pt", "subleading lepton pt", 50, 0, 200);
00135     subsubleading_l_pt = dbe->book1D("subsubleading_l_pt", "subsubleading lepton pt", 50, 0, 200);
00136     leading_l_eta = dbe->book1D("leading_l_eta", "leading lepton eta",  30, -3.5, 3.5);
00137     subleading_l_eta = dbe->book1D("subleading_l_eta", "subleading lepton eta",  30, -3.5, 3.5);
00138     subsubleading_l_eta = dbe->book1D("subsubleading_l_eta", "subsubleading lepton eta",  30, -3.5, 3.5);
00139     mll= dbe->book1D("mll", "ll mass (all combinations)", 50, 0, 200);
00140     ptll= dbe->book1D("ptll", "ll Transverse Momentum (all combinations)", 50, 0, 200);
00141     mlll= dbe->book1D("mlll", "lll mass ", 50, 0, 200);
00142     ptlll= dbe->book1D("ptlll", "lll Transverse Momentum ", 50, 0, 2000);
00143     mlllnununu= dbe->book1D("mlllnununu", "lll nununu mass ", 50, 0, 2000);
00144     mtlllnununu= dbe->book1D("mtlllnununu", "lll nununu transverse mass ", 50, 0, 2000);
00145     ptlllnununu= dbe->book1D("ptlllnununu", "lll nununu Transverse Momentum ", 50, 0, 2000);
00146           
00147 
00148   }
00149   return;
00150 }
00151 
00152 void VVVValidation::endJob(){return;}
00153 void VVVValidation::beginRun(const edm::Run& iRun,const edm::EventSetup& iSetup)
00154 {
00155   iSetup.getData( fPDGTable );
00156   return;
00157 }
00158 void VVVValidation::endRun(const edm::Run& iRun,const edm::EventSetup& iSetup){return;}
00159 void VVVValidation::analyze(const edm::Event& iEvent,const edm::EventSetup& iSetup)
00160 { 
00161 
00162   unsigned int initSize = 1000;
00163 
00164   edm::Handle<HepMCProduct> evt;
00165   iEvent.getByLabel(hepmcCollection_, evt);
00166 
00167   HepMC::GenEvent *myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));
00168 
00169   double weight = _wmanager.weight(iEvent);
00170 
00171   nEvt->Fill(0.5, weight);
00172 
00173   int st_3_leptons=0;
00174   int n_particles=0;
00175   int Wpst3=0;
00176   int Wmst3=0;
00177   int taust3=0;
00178 
00179   std::vector<const HepMC::GenParticle*> mothercoll;
00180   std::vector<const HepMC::GenParticle*> GenLeptons;
00181   std::vector<const HepMC::GenParticle*> GenNeutrinos;
00182   mothercoll.reserve(initSize);
00183   GenLeptons.reserve(initSize);
00184   GenNeutrinos.reserve(initSize);
00185   for (HepMC::GenEvent::particle_const_iterator iter = myGenEvent->particles_begin(); iter != myGenEvent->particles_end(); ++iter){
00186     double id = (*iter)->pdg_id();
00187 
00188     if((*iter)->status()==3&&(fabs(id)==11 ||fabs(id)==13) && (*iter)->momentum().perp()>20.&& fabs((*iter)->momentum().eta())<2.5){
00189       st_3_leptons++;
00190     }
00191     if((*iter)->status()==1&&(fabs(id)==16))taust3++;
00192     if((*iter)->status()==3&&(id==24))Wpst3++;
00193     if((*iter)->status()==3&&(id==-24))Wmst3++;
00194     if((*iter)->status()==3&&(fabs(id)==24 || fabs(id)==23)){
00195      mothercoll.push_back(*iter);
00196     }
00197 
00198     if((*iter)->status()==_lepStatus && (*iter)->momentum().perp()>20.&& fabs((*iter)->momentum().eta())<2.5 ){
00199       n_particles ++;
00200       if (fabs(id)==11 ||fabs(id)==13){
00201         GenLeptons.push_back(*iter);
00202       }
00203       if (fabs(id)==12 ||fabs(id)==14){
00204         GenNeutrinos.push_back(*iter);
00205       }
00206     }
00207   }
00208 
00209   vector<TLorentzVector> wpluss;
00210   vector<TLorentzVector> wmins;
00211   vector<TLorentzVector> z;
00212   wpluss.clear();
00213   wmins.clear();
00214   z.clear();
00215   int Wmin=0;
00216   int Wpl=0;
00217   int Z=0;
00218   int l_bcode[3];
00219   int nu_bcode[3];
00220 
00221 
00222   if(taust3>0)return;
00223   
00224   for(unsigned int i = 0 ; i< mothercoll.size();i++){
00225     double id = mothercoll[i]->pdg_id();
00226 
00228     //&&&&W bosons&&&&//
00230     if(fabs(id)==24){
00231       double dr_min=999.;
00232       double dr=999.;
00233       for(unsigned int k=0 ; k< GenLeptons.size() ;k++){
00234         for(unsigned int kl=0 ; kl< GenNeutrinos.size() ;kl++){
00235           double lepton_px=GenLeptons[k]->momentum().px();
00236           double lepton_py=GenLeptons[k]->momentum().py();
00237           double lepton_pz=GenLeptons[k]->momentum().pz();
00238           double lepton_e=GenLeptons[k]->momentum().e();
00239           TLorentzVector l;
00240           l.SetPxPyPzE(lepton_px,lepton_py,lepton_pz,lepton_e); 
00241           double nu_px=GenNeutrinos[kl]->momentum().px();
00242           double nu_py=GenNeutrinos[kl]->momentum().py();
00243           double nu_pz=GenNeutrinos[kl]->momentum().pz();
00244           double nu_e=GenNeutrinos[kl]->momentum().e();
00245           TLorentzVector nu;
00246           nu.SetPxPyPzE(nu_px,nu_py,nu_pz,nu_e);
00247           double l_id= GenLeptons[k]->pdg_id();  
00248           double nu_id= GenNeutrinos[kl]->pdg_id();
00249           dr= deltaR((l+nu).PseudoRapidity(),(l+nu).Phi(),mothercoll[i]->momentum().eta(),mothercoll[i]->momentum().phi());
00250           if((id*l_id)<0 &&(l_id*nu_id)<0 &&( fabs(nu_id)== (fabs(l_id)+1) )&& dr<0.5/*&& (l+nu).M()>6.*/){
00251             if(dr<dr_min)dr_min=dr;
00252             if(dr>dr_min)continue;
00253           }
00254         }
00255       }
00256       for(unsigned int k=0 ; k< GenLeptons.size() ;k++){
00257         for(unsigned int kl=0 ; kl< GenNeutrinos.size() ;kl++){
00258           double lepton_px=GenLeptons[k]->momentum().px();
00259           double lepton_py=GenLeptons[k]->momentum().py();
00260           double lepton_pz=GenLeptons[k]->momentum().pz();
00261           double lepton_e=GenLeptons[k]->momentum().e();
00262           TLorentzVector l;
00263           l.SetPxPyPzE(lepton_px,lepton_py,lepton_pz,lepton_e); 
00264           double nu_px=GenNeutrinos[kl]->momentum().px();
00265           double nu_py=GenNeutrinos[kl]->momentum().py();
00266           double nu_pz=GenNeutrinos[kl]->momentum().pz();
00267           double nu_e=GenNeutrinos[kl]->momentum().e();
00268           TLorentzVector nu;
00269           nu.SetPxPyPzE(nu_px,nu_py,nu_pz,nu_e);
00270           double l_id= GenLeptons[k]->pdg_id();  
00271           double nu_id= GenNeutrinos[kl]->pdg_id();
00272           double der= deltaR((l+nu).PseudoRapidity(),(l+nu).Phi(),mothercoll[i]->momentum().eta(),mothercoll[i]->momentum().phi());
00273           if((id*l_id)<0 && (l_id*nu_id)<0 &&( fabs(nu_id)== (fabs(l_id)+1) )&& der==dr_min){
00274             l_bcode[i]=GenLeptons[k]->barcode();
00275             nu_bcode[i]=GenNeutrinos[kl]->barcode();
00276             if((i==0)|| (i==1 &&  (l_bcode[i]!=l_bcode[i-1]) && (nu_bcode[i]!=nu_bcode[i-1])  )||
00277                (i==2 &&  (l_bcode[i]!=l_bcode[i-1]) && (nu_bcode[i]!=nu_bcode[i-1]) && (l_bcode[i]!=l_bcode[i-2]) && (nu_bcode[i]!=nu_bcode[i-2]) )
00278               ){
00279               if(id==24){
00280                 Wpl++;
00281                 wpluss.push_back((l+nu));
00282                 h_mWplus->Fill((l+nu).M(),weight);
00283                 h_phiWplus->Fill((l+nu).Phi(),weight);
00284                 h_ptWplus->Fill((l+nu).Pt(),weight);
00285                 h_yWplus->Fill((l+nu).Rapidity(),weight);
00286               }
00287               if(id==-24){
00288                 Wmin++;
00289                 wmins.push_back((l+nu));
00290                 h_mWminus->Fill((l+nu).M(),weight);
00291                 h_phiWminus->Fill((l+nu).Phi(),weight);
00292                 h_ptWminus->Fill((l+nu).Pt(),weight);
00293                 h_yWminus->Fill((l+nu).Rapidity(),weight);
00294               }
00295             } 
00296           }
00297         }
00298       }
00299     }
00300     
00302     //&&&&Z bosons&&&&//
00304     if(fabs(id)==23){
00305       double dr_min=999.;
00306       double dr=999.;
00307       for(unsigned int k=0 ; k< GenLeptons.size() ;k++){
00308         for(unsigned int kl=k ; kl< GenLeptons.size() ;kl++){
00309           if(k==kl)continue;
00310           double lepton_px=GenLeptons[k]->momentum().px();
00311           double lepton_py=GenLeptons[k]->momentum().py();
00312           double lepton_pz=GenLeptons[k]->momentum().pz();
00313           double lepton_e=GenLeptons[k]->momentum().e();
00314           TLorentzVector l;
00315           l.SetPxPyPzE(lepton_px,lepton_py,lepton_pz,lepton_e); 
00316           double nu_px=GenLeptons[kl]->momentum().px();
00317           double nu_py=GenLeptons[kl]->momentum().py();
00318           double nu_pz=GenLeptons[kl]->momentum().pz();
00319           double nu_e=GenLeptons[kl]->momentum().e();
00320           TLorentzVector nu;
00321           nu.SetPxPyPzE(nu_px,nu_py,nu_pz,nu_e);
00322           double l_id= GenLeptons[k]->pdg_id();  
00323           double nu_id= GenLeptons[kl]->pdg_id();
00324           dr= deltaR((l+nu).PseudoRapidity(),(l+nu).Phi(),mothercoll[i]->momentum().eta(),mothercoll[i]->momentum().phi());
00325           if((l_id*nu_id)<0 &&( fabs(nu_id)== (fabs(l_id)) && dr<0.5)/*&& (l+nu).M()>6.*/){
00326             if(dr<dr_min)dr_min=dr;
00327             if(dr>dr_min)continue;
00328           }
00329         }
00330       }
00331       for(unsigned int k=0 ; k< GenLeptons.size() ;k++){
00332         for(unsigned int kl=k ; kl< GenLeptons.size() ;kl++){
00333           if(k==kl)continue;
00334           double lepton_px=GenLeptons[k]->momentum().px();
00335           double lepton_py=GenLeptons[k]->momentum().py();
00336           double lepton_pz=GenLeptons[k]->momentum().pz();
00337           double lepton_e=GenLeptons[k]->momentum().e();
00338           TLorentzVector l;
00339           l.SetPxPyPzE(lepton_px,lepton_py,lepton_pz,lepton_e); 
00340           double nu_px=GenLeptons[kl]->momentum().px();
00341           double nu_py=GenLeptons[kl]->momentum().py();
00342           double nu_pz=GenLeptons[kl]->momentum().pz();
00343           double nu_e=GenLeptons[kl]->momentum().e();
00344           TLorentzVector nu;
00345           nu.SetPxPyPzE(nu_px,nu_py,nu_pz,nu_e);
00346           double l_id= GenLeptons[k]->pdg_id();  
00347           double nu_id= GenLeptons[kl]->pdg_id();
00348           double der= deltaR((l+nu).PseudoRapidity(),(l+nu).Phi(),mothercoll[i]->momentum().eta(),mothercoll[i]->momentum().phi());
00349           if((l_id*nu_id)<0 &&( fabs(nu_id)== (fabs(l_id)) )&& der==dr_min ){
00350             l_bcode[i]=GenLeptons[k]->barcode();
00351             nu_bcode[i]=GenLeptons[kl]->barcode();
00352            if((i==0)|| (i==1 &&  (l_bcode[i]!=l_bcode[i-1]) && (nu_bcode[i]!=nu_bcode[i-1])  )||
00353                (i==2 &&  (l_bcode[i]!=l_bcode[i-1]) && (nu_bcode[i]!=nu_bcode[i-1]) && (l_bcode[i]!=l_bcode[i-2]) && (nu_bcode[i]!=nu_bcode[i-2]) )
00354              ){
00355               Z++;
00356               z.push_back((l+nu));
00357               h_mZ->Fill((l+nu).M(),weight);
00358               h_phiZ->Fill((l+nu).Phi(),weight);
00359               h_ptZ->Fill((l+nu).Pt(),weight);
00360               h_yZ->Fill((l+nu).Rapidity(),weight); 
00361             }          
00362           }
00363         }
00364       }
00365     }
00366   }
00367 
00368 
00369   if((Wmin+Wpl)>3) cout<<"3ten fazla W adayı?!?"<<endl;
00370 
00371   if( ((Wmin+Wpl)==3) || ((Wmin+Wpl)==2 && Z==1) || (( Wmin+Wpl )==1 && Z==2) || (Z==3) ){
00372 
00373 
00374     for(unsigned int i=0; i<wmins.size();i++){
00375       h_mWminus_3b->Fill(wmins[i].M(),weight);
00376       h_phiWminus_3b->Fill(wmins[i].Phi(),weight);
00377       h_ptWminus_3b->Fill(wmins[i].Pt(),weight);
00378       h_yWminus_3b->Fill(wmins[i].Rapidity(),weight);
00379     }
00380     for(unsigned int i=0; i<wpluss.size();i++){
00381       h_mWplus_3b->Fill(wpluss[i].M(),weight);
00382       h_phiWplus_3b->Fill(wpluss[i].Phi(),weight);
00383       h_ptWplus_3b->Fill(wpluss[i].Pt(),weight);
00384       h_yWplus_3b->Fill(wpluss[i].Rapidity(),weight);
00385     }
00386     for(unsigned int i=0; i<z.size();i++){
00387       h_mZ_3b->Fill(z[i].M(),weight);
00388       h_phiZ_3b->Fill(z[i].Phi(),weight);
00389       h_ptZ_3b->Fill(z[i].Pt(),weight);
00390       h_yZ_3b->Fill(z[i].Rapidity(),weight);
00391     }
00392  
00393     if((Wmin+Wpl)==3){
00394         TLorentzVector W1;
00395         TLorentzVector W2;
00396         TLorentzVector W3;
00397       if(Wmin==2 && Wpl==1){
00398         W1=wmins[0];
00399         W2=wmins[1];
00400         W3=wpluss[0];
00401       }
00402       if(Wmin==1 && Wpl==2){
00403         W1=wmins[0];
00404         W2=wpluss[0];
00405         W3=wpluss[1];
00406       }
00407       if(W1.M()<10. ||W2.M()<10.||W3.M()<10.)cout<<taust3<<endl;
00408       h_mWWW->Fill((W1+W2+W3).M(),weight);
00409       h_phiWWW->Fill((W1+W2+W3).Phi(),weight);
00410       h_ptWWW->Fill((W1+W2+W3).Pt(),weight);
00411       h_yWWW->Fill((W1+W2+W3).Rapidity(),weight);
00412 
00413       h_mWW->Fill((W1+W2).M(),weight);
00414       h_mWW->Fill((W1+W3).M(),weight);
00415       h_mWW->Fill((W2+W3).M(),weight);
00416 
00417       h_phiWW->Fill((W1+W2).Phi(),weight);
00418       h_phiWW->Fill((W1+W3).Phi(),weight);
00419       h_phiWW->Fill((W2+W3).Phi(),weight);
00420 
00421       h_ptWW->Fill((W1+W2).Pt(),weight);
00422       h_ptWW->Fill((W1+W3).Pt(),weight);
00423       h_ptWW->Fill((W2+W3).Pt(),weight);
00424 
00425       h_yWW->Fill((W1+W2).Rapidity(),weight);
00426       h_yWW->Fill((W1+W3).Rapidity(),weight);
00427       h_yWW->Fill((W2+W3).Rapidity(),weight);
00428 
00429 
00430     }
00431     if((Wmin+Wpl)==2 && Z==1){
00432       TLorentzVector W1;
00433       TLorentzVector W2;
00434       TLorentzVector Z1;
00435       Z1=z[0];
00436       if(Wmin==1 &&Wpl==1){
00437         W1=wmins[0];
00438         W2=wpluss[0];
00439       }  
00440       if(Wmin==2 &&Wpl==0){
00441         W1=wmins[0];
00442         W2=wmins[1];
00443       }  
00444       if(Wmin==0 &&Wpl==2){
00445         W1=wpluss[0];
00446         W2=wpluss[1];
00447       }
00448       h_mWW->Fill((W1+W2).M(),weight);
00449       h_phiWW->Fill((W1+W2).Phi(),weight);
00450       h_ptWW->Fill((W1+W2).Pt(),weight);
00451       h_yWW->Fill((W1+W2).Rapidity(),weight);
00452 
00453       h_mWZ->Fill((W1+Z1).M(),weight);
00454       h_mWZ->Fill((W2+Z1).M(),weight);
00455 
00456       h_phiWZ->Fill((W1+Z1).Phi(),weight);
00457       h_phiWZ->Fill((W2+Z1).Phi(),weight);
00458 
00459       h_ptWZ->Fill((W1+Z1).Pt(),weight);
00460       h_ptWZ->Fill((W2+Z1).Pt(),weight);
00461 
00462       h_yWZ->Fill((W1+Z1).Rapidity(),weight);
00463       h_yWZ->Fill((W2+Z1).Rapidity(),weight);
00464 
00465       h_mWWZ->Fill((W1+W2+Z1).M(),weight);
00466       h_phiWWZ->Fill((W1+W2+Z1).Phi(),weight);
00467 
00468       h_ptWWZ->Fill((W1+W2+Z1).Pt(),weight);
00469       h_yWWZ->Fill((W1+W2+Z1).Rapidity(),weight);
00470 
00471                   
00472     }
00473 
00474     if((Wmin+Wpl)==1 && Z==2){
00475       TLorentzVector Z1;
00476       TLorentzVector Z2;
00477       TLorentzVector W1;
00478       Z1=z[0];
00479       Z2=z[1];
00480       if(Wmin==1){
00481         W1=wmins[0];
00482       }  
00483       if(Wpl==1){
00484         W1=wpluss[0];
00485       }  
00486       h_mZZ->Fill((Z1+Z2).M(),weight);
00487       h_phiZZ->Fill((Z1+Z2).Phi(),weight);
00488       h_ptZZ->Fill((Z1+Z2).Pt(),weight);
00489       h_yZZ->Fill((Z1+Z2).Rapidity(),weight);
00490 
00491       h_mWZ->Fill((W1+Z1).M(),weight);
00492       h_mWZ->Fill((W1+Z2).M(),weight);
00493 
00494       h_phiWZ->Fill((W1+Z1).Phi(),weight);
00495       h_phiWZ->Fill((W1+Z2).Phi(),weight);
00496 
00497       h_ptWZ->Fill((W1+Z1).Pt(),weight);
00498       h_ptWZ->Fill((W1+Z2).Pt(),weight);
00499 
00500       h_yWZ->Fill((W1+Z1).Rapidity(),weight);
00501       h_yWZ->Fill((W1+Z2).Rapidity(),weight);
00502 
00503       h_mWZZ->Fill((Z1+Z2+W1).M(),weight);
00504       h_phiWZZ->Fill((Z1+Z2+W1).Phi(),weight);
00505       h_ptWZZ->Fill((Z1+Z2+W1).Pt(),weight);
00506       h_yWZZ->Fill((Z1+Z2+W1).Rapidity(),weight);
00507 
00508                   
00509     }
00510 
00511     if(Z==3){
00512       TLorentzVector Z1;
00513       TLorentzVector Z2;
00514       TLorentzVector Z3;
00515 
00516       Z1=z[0];
00517       Z2=z[1];
00518       Z3=z[2];
00519 
00520       if(Z1.M()<10. ||Z2.M()<10.||Z3.M()<10.)cout<<taust3<<endl;
00521 
00522       h_mZZZ->Fill((Z1+Z2+Z3).M(),weight);
00523       h_phiZZZ->Fill((Z1+Z2+Z3).Phi(),weight);
00524       h_ptZZZ->Fill((Z1+Z2+Z3).Pt(),weight);
00525       h_yZZZ->Fill((Z1+Z2+Z3).Rapidity(),weight);
00526 
00527       h_mZZ->Fill((Z1+Z2).M(),weight);
00528       h_mZZ->Fill((Z1+Z3).M(),weight);
00529       h_mZZ->Fill((Z2+Z3).M(),weight);
00530 
00531       h_phiZZ->Fill((Z1+Z2).Phi(),weight);
00532       h_phiZZ->Fill((Z1+Z3).Phi(),weight);
00533       h_phiZZ->Fill((Z2+Z3).Phi(),weight);
00534 
00535       h_ptZZ->Fill((Z1+Z2).Pt(),weight);
00536       h_ptZZ->Fill((Z1+Z3).Pt(),weight);
00537       h_ptZZ->Fill((Z2+Z3).Pt(),weight);
00538 
00539       h_yZZ->Fill((Z1+Z2).Rapidity(),weight);
00540       h_yZZ->Fill((Z1+Z3).Rapidity(),weight);
00541       h_yZZ->Fill((Z2+Z3).Rapidity(),weight);
00542 
00543     }
00544 
00545 
00546           
00547     edm::Handle<reco::GenJetCollection> genJets;
00548     iEvent.getByLabel(genjetCollection_, genJets );
00549 
00550     std::vector<const reco::GenJet*> selected_jets;
00551     selected_jets.reserve(initSize);
00552 
00553     int nJets = 0;
00554     int nJetso1 = 0;
00555     int nJetso30 = 0;
00556     int nJetso50 = 0;
00557     int nJetso100 = 0;
00558     int nJetsCentral = 0;
00559     double totPt = 0.;
00560     double dr=99.;
00561     std::vector<double> jetEta;
00562     jetEta.reserve(initSize);
00563     for (reco::GenJetCollection::const_iterator iter=genJets->begin();iter!=genJets->end();++iter){
00564       dr=0;
00565       bool matched_lepton=false;
00566       for(unsigned int i=0 ; i< GenLeptons.size() ;i++){
00567         dr= deltaR(GenLeptons[i]->momentum().eta(),GenLeptons[i]->momentum().phi(),(*iter).eta(),(*iter).phi());
00568         h_dr->Fill(dr,weight);
00569         if(dr<0.5)matched_lepton=true;
00570       }
00571       if(matched_lepton)continue;
00572       selected_jets.push_back(&*iter);
00573       nJets++;
00574       double pt = (*iter).pt();
00575       totPt += pt;
00576       if (pt > 1.) nJetso1++;
00577       double eta = (*iter).eta();
00578       if (std::fabs(eta) < 5.&&pt > 30.) nJetso30++;
00579       if (std::fabs(eta) < 5.&&pt > 50.) nJetso50++;
00580       if (std::fabs(eta) < 5.&&pt > 100.) nJetso100++;
00581       if ( std::fabs(eta) < 2.5 ) nJetsCentral++;
00582       jetEta.push_back(eta);
00583 
00584       genJetEnergy->Fill(std::log10((*iter).energy()),weight);
00585       genJetPt->Fill(std::log10(pt),weight);
00586       genJetEta->Fill(eta,weight);
00587       genJetPhi->Fill((*iter).phi()/CLHEP::degree,weight);
00588     }
00589     if(nJetso30==2){
00590  
00591       TLorentzVector j1;
00592       TLorentzVector j2;
00593       j1.SetPtEtaPhiE(selected_jets[0]->pt(),selected_jets[0]->eta(),selected_jets[0]->phi(),selected_jets[0]->energy());
00594       j2.SetPtEtaPhiE(selected_jets[1]->pt(),selected_jets[1]->eta(),selected_jets[1]->phi(),selected_jets[1]->energy());
00595       WW_TwoJEt_JetsM->Fill((j1+j2).M(),weight);
00596     } 
00597     if(nJetso30>0)h_l_jet_eta->Fill(selected_jets[0]->eta());
00598     if(nJetso30>0)h_l_jet_pt->Fill(selected_jets[0]->pt());
00599     if(nJetso30>1)h_sl_jet_eta->Fill(selected_jets[1]->eta());
00600     if(nJetso30>1)h_sl_jet_pt->Fill(selected_jets[1]->pt());
00601     if(nJetso30>2)h_ssl_jet_eta->Fill(selected_jets[2]->eta());
00602     if(nJetso30>2)h_ssl_jet_pt->Fill(selected_jets[2]->pt());
00603 
00604     genJetMult->Fill(nJets,weight);
00605     genJetPto1->Fill(nJetso1,weight);
00606     genJetPto30->Fill(nJetso30,weight);
00607     genJetPto50->Fill(nJetso50,weight);
00608     genJetPto100->Fill(nJetso100,weight);
00609     genJetCentral->Fill(nJetsCentral,weight);
00610 
00611     genJetTotPt->Fill(totPt,weight);
00612 
00613     double deltaEta = 999.;
00614     if ( jetEta.size() > 1 ) {
00615       for (unsigned int i = 0; i < jetEta.size(); i++){
00616         for (unsigned int j = i+1; j < jetEta.size(); j++){
00617           deltaEta = std::min(deltaEta,std::fabs(jetEta[i]-jetEta[j]));
00618         }
00619       }
00620     }
00621 
00622     genJetDeltaEtaMin->Fill(deltaEta,weight);
00623   }
00624 
00625   if(GenLeptons.size()>0 && GenNeutrinos.size()>0 ){
00626     std::sort(GenLeptons.begin(), GenLeptons.end(), HepMCValidationHelper::sortByPt); 
00627     std::sort(GenNeutrinos.begin(), GenNeutrinos.end(), HepMCValidationHelper::sortByPt);
00628     if(GenLeptons.size()>0)    leading_l_pt->Fill(GenLeptons[0]->momentum().perp(),weight);
00629     if(GenLeptons.size()>1)    subleading_l_pt->Fill(GenLeptons[1]->momentum().perp(),weight);
00630     if(GenLeptons.size()>2)    subsubleading_l_pt->Fill(GenLeptons[2]->momentum().perp(),weight);
00631     if(GenLeptons.size()>0)    leading_l_eta->Fill(GenLeptons[0]->momentum().eta(),weight);
00632     if(GenLeptons.size()>1)    subleading_l_eta->Fill(GenLeptons[1]->momentum().eta(),weight);
00633     if(GenLeptons.size()>2)    subsubleading_l_eta->Fill(GenLeptons[2]->momentum().eta(),weight);
00634   }
00635   if(GenLeptons.size()>1 ){
00636     for(unsigned int i = 0; i<GenLeptons.size();i++){
00637       for(unsigned int j = i; j<GenLeptons.size();j++){
00638         if(j==i)continue;
00639         TLorentzVector lep1(GenLeptons[i]->momentum().x(), GenLeptons[i]->momentum().y(), GenLeptons[i]->momentum().z(), GenLeptons[i]->momentum().t()); 
00640         TLorentzVector lep2(GenLeptons[j]->momentum().x(), GenLeptons[j]->momentum().y(), GenLeptons[j]->momentum().z(), GenLeptons[j]->momentum().t()); 
00641         mll->Fill((lep1+lep2).M(),weight);
00642         ptll->Fill((lep1+lep2).Pt(),weight);
00643       }
00644     }
00645   }
00646   if(GenLeptons.size()>2 && GenNeutrinos.size()>2 ){
00647     TLorentzVector lep1(GenLeptons[0]->momentum().x(), GenLeptons[0]->momentum().y(), GenLeptons[0]->momentum().z(), GenLeptons[0]->momentum().t()); 
00648     TLorentzVector lep2(GenLeptons[1]->momentum().x(), GenLeptons[1]->momentum().y(), GenLeptons[1]->momentum().z(), GenLeptons[1]->momentum().t()); 
00649     TLorentzVector lep3(GenLeptons[2]->momentum().x(), GenLeptons[2]->momentum().y(), GenLeptons[2]->momentum().z(), GenLeptons[2]->momentum().t());
00650     TLorentzVector nu1(GenNeutrinos[0]->momentum().x(), GenNeutrinos[0]->momentum().y(), GenNeutrinos[0]->momentum().z(), GenNeutrinos[0]->momentum().t()); 
00651     TLorentzVector nu2(GenNeutrinos[1]->momentum().x(), GenNeutrinos[1]->momentum().y(), GenNeutrinos[1]->momentum().z(), GenNeutrinos[1]->momentum().t()); 
00652     TLorentzVector nu3(GenNeutrinos[2]->momentum().x(), GenNeutrinos[2]->momentum().y(), GenNeutrinos[2]->momentum().z(), GenNeutrinos[2]->momentum().t()); 
00653     mlll->Fill((lep1+lep2+lep3).M(),weight);
00654     ptlll->Fill((lep1+lep2+lep3).Pt(),weight);
00655     mlllnununu->Fill((lep1+lep2+lep3+nu1+nu2+nu3).M(),weight);
00656     mtlllnununu->Fill((lep1+lep2+lep3+nu1+nu2+nu3).Mt(),weight);
00657     ptlllnununu->Fill((lep1+lep2+lep3+nu1+nu2+nu3).Pt(),weight);
00658                 
00659   }
00660         
00661   delete myGenEvent;
00662 }
00663 int VVVValidation::getParentBarcode(HepMC::GenParticle* it)
00664 {
00665     int id = 0;
00666     if ( (it)->production_vertex() && (it)-> status()==3) {
00667         for ( HepMC::GenVertex::particle_iterator mother 
00668                   = (it)->production_vertex()->particles_begin(HepMC::parents);mother != (it)->production_vertex()->particles_end(HepMC::parents); ++mother ) {
00669 
00670            if((fabs((*mother)->pdg_id())==24)) id = (*mother)->barcode();
00671         }
00672     }
00673     return id;
00674 }
00675 #include "FWCore/PluginManager/interface/ModuleDef.h"
00676 #include "FWCore/Framework/interface/MakerMacros.h"
00677 
00678 DEFINE_FWK_MODULE (VVVValidation);
00679
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/Validation/EventGenerator/plugins/VVVValidation.cc
