/*======================== Au Diya 27 October 2020 *======================== // fighting corona // **************************** Analysis strategy is to get the hist for 4mu 4e 2mu2e ****************************/ // system include files #include #include #include #include #include #include #include #include "TSystem.h" #include "math.h" // for Root histogram #include "TH1.h" #include "TTree.h" #include "TH2.h" #include "TLorentzVector.h" // for Delphes #ifdef __CLING__ R__LOAD_LIBRARY(libDelphes) #include "classes/DelphesClasses.h" #include "external/ExRootAnalysis/ExRootTreeReader.h" #include "external/ExRootAnalysis/ExRootResult.h" #else class ExRootTreeReader; class ExRootResult; #endif template void CollectionFilter(const TClonesArray& inColl ,vector& outColl, Double_t ptMin=30, Double_t etaMax=2.5) { const TObject *object; for (Int_t i = 0; i < inColl.GetEntriesFast(); i++) { object = inColl.At(i); const T *t = static_cast(object); if(t->P4().Pt() < ptMin) continue; if(TMath::Abs(t->P4().Eta()) > etaMax) continue; outColl.push_back(t); } } //------------------------------------------------------------------------ void Audiya4lAnalyzer(const char *fileName="/Users/applestudio/Desktop/sample28sept/ALLROOTFILES"){ gSystem->Load("libDelphes"); // Create chain of root trees TChain chain("Delphes"); //4e chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/e10.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/e15.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/e20.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/e25.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/e30.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/e35.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/e40.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/e45.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/e50.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/e55.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/e60.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/e65.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/e70.root"); //4mu chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/mu10.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/mu15.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/mu20.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/mu25.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/mu30.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/mu35.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/mu40.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/mu45.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/mu50.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/mu55.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/mu60.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/mu65.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/mu70.root"); //2mu2e chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2mu2e10.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2mu2e15.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2mu2e20.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2mu2e25.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2mu2e30.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2mu2e35.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2mu2e40.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2mu2e45.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2mu2e50.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2mu2e55.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2mu2e60.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2mu2e65.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2mu2e70.root"); //2e2mu chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2e2mu10.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2e2mu15.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2e2mu20.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2e2mu25.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2e2mu30.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2e2mu35.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2e2mu40.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2e2mu45.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2e2mu50.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2e2mu55.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2e2mu60.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2e2mu65.root"); chain.Add("/Users/applestudio/Desktop/sample28sept/ALLROOTFILES/2e2mu70.root"); // Create object of class ExRootTreeReader ExRootTreeReader *treeReader = new ExRootTreeReader(&chain); Long64_t numberOfEntries = treeReader->GetEntries(); // Get pointers to branches used in this analysis TClonesArray *branchParticle = treeReader->UseBranch("Particle"); TClonesArray *branchElectron = treeReader->UseBranch("Electron"); TClonesArray *branchMuon = treeReader->UseBranch("Muon"); // Declare variable const Electron *elec1, *elec2, *elec3, *elec4; const Muon *muon1, *muon2, *muon3, *muon4; vector *electrons = new vector();; vector *muons = new vector();; //TLorentzVector ele; //TLorentzVector mu; //ele.SetPtEtaPhiM(ele->PT, ele->Eta, ele->Phi, mass_ele); //mu.SetPtEtaPhiM(mu->PT, mu->Eta, mu->Phi, mass_mu); //MuonVector muons(PT, Eta, Phi, Charge, E); //ElectronVector electrons(PT, Eta, Phi, Charge, E); Double_t mass, massEE, massMM; double mZ12, mZ34, mZ13, mZ24, mZ14, mZ23; double mass4mu; double mZa, mZb; double mZ, eZ12, eZ34, eZ13, eZ24, eZ14, eZ23; double pxZ12, pxZ34, pxZ13, pxZ24, pxZ14, pxZ23; double pyZ12, pyZ34, pyZ13, pyZ24, pyZ14, pyZ23; double pzZ12, pzZ34, pzZ13, pzZ24, pzZ14, pzZ23; double pZ12, pZ34, pZ13, pZ24, pZ14, pZ23; double pTZ12, pTZ34, pTZ13, pTZ24, pTZ14, pTZ23; double dZ12, dZ34, dZ13, dZ24, dZ14, dZ23; double dZc1, dZc2, dZc3; double eZa, pxZa, pyZa, pzZa, pTZa; double eZb, pxZb, pyZb, pzZb, pTZb; double mass4e; double px_e1, px_e2, px_e3, px_e4; double py_e1, py_e2, py_e3, py_e4; double pz_e1, pz_e2, pz_e3, pz_e4; double E_e1, E_e2, E_e3, E_e4; double mass2mu2e, pt_2mu2e, eta_2mu2e, phi_2mu2e; double px2mu2e, py2mu2e, pz2mu2e, E2mu2e; double pt_2mu1, pt_2mu2, pt_2e1, pt_2e2; double eta_2mu1, eta_2mu2, eta_2e1, eta_2e2; double phi_2mu1, phi_2mu2, phi_2e1, phi_2e2; int cas_2mu1, cas_2mu2, cas_2e1, cas_2e2; double px_2mu1, px_2mu2, px_2e1, px_2e2; double py_2mu1, py_2mu2, py_2e1, py_2e2; double pz_2mu1, pz_2mu2, pz_2e1, pz_2e2; double E_2mu1, E_2mu2, E_2e1, E_2e2; TLorentzVector p4Za, p4Zb, p4H; // Book histograms & set axis lable //**Dimuon // ZTo2mu mass TH1D *h_mZ_2mu = new TH1D("massZto2muon", "mass of Z to 2 muon",120, 40., 120.); h_mZ_2mu->GetXaxis()->SetTitle("Invariant mass of dimuon (in GeV/c^2)"); h_mZ_2mu->GetYaxis()->SetTitle("Number of Events"); // ZTo2e mass TH1D *h_mZ_2e = new TH1D("massZto2e", "mass of Z to 2e", 120, 40., 120.); h_mZ_2e->GetXaxis()->SetTitle("Invariant Mass for Nelectron=2 (in GeV/c^2)"); h_mZ_2e->GetYaxis()->SetTitle("Number of Events"); // These histograms are for 4 MUON reconstruction with different combinations // First combination: 1234 // Mass Z12 TH1D *h_mZ12_4mu = new TH1D("mZ12_4mu", "mass of Z12", 75, 0., 150.); h_mZ12_4mu->GetXaxis()->SetTitle("Invariant mass of dimuon (in GeV/c^2)"); h_mZ12_4mu->GetYaxis()->SetTitle("Number of Events"); // Mass Z34 TH1D *h_mZ34_4mu = new TH1D("mZ34_4mu", "mass of Z34", 75, 0., 150.); h_mZ34_4mu->GetXaxis()->SetTitle("Invariant mass of dimuon (in GeV/c^2)"); h_mZ34_4mu->GetYaxis()->SetTitle("Number of Events"); // Second combination: 1324 // Mass Z13 TH1D *h_mZ13_4mu = new TH1D("mZ13_4mu", "mass of Z13", 75, 0., 150.); h_mZ13_4mu->GetXaxis()->SetTitle("Invariant mass of dimuon (in GeV/c^2)"); h_mZ13_4mu->GetYaxis()->SetTitle("Number of Events"); // Mass Z24 TH1D *h_mZ24_4mu = new TH1D("mZ24_4mu", "mass of Z24", 75, 0., 150.); h_mZ24_4mu->GetXaxis()->SetTitle("Invariant mass of dimuon"); h_mZ24_4mu->GetYaxis()->SetTitle("Number of Events"); // Third combination: 1423 // Mass Z14 TH1D *h_mZ14_4mu = new TH1D("mZ14_4mu", "mass of Z14", 75, 0., 150.); h_mZ14_4mu->GetXaxis()->SetTitle("Invariant mass of dimuon (in GeV/c^2)"); h_mZ14_4mu->GetYaxis()->SetTitle("Number of Events"); // Mass Z23 TH1D *h_mZ23_4mu = new TH1D("mZ23_4mu", "mass of Z23", 75, 0., 150.); h_mZ23_4mu->GetXaxis()->SetTitle("Invariant mass of dimuon (in GeV/c^2)"); h_mZ23_4mu->GetYaxis()->SetTitle("Number of Events"); // Mass Za: mass of ZTo2mu closest to Z mass TH1D *h_mZa_4mu = new TH1D("mZa_4mu", "mass Za", 120, 0., 120.); h_mZa_4mu->GetXaxis()->SetTitle("Invariant mass of dimuon (in GeV/c^2)"); h_mZa_4mu->GetYaxis()->SetTitle("Number of Events"); // Mass Zb: mass of ZTo2mu not closest to Z mass TH1D *h_mZb_4mu = new TH1D("mZb_4mu", "mass Zb", 120, 0., 120.); h_mZb_4mu->GetXaxis()->SetTitle("Invariant mass of dimuon (in GeV/c^2)"); h_mZb_4mu->GetYaxis()->SetTitle("Number of Events"); // 4muon mass o muon (full mass range) TH1D *h_m4_m4mu = new TH1D("mass4mu_full", "mass of 4 muon", 300, 0., 900.); h_m4_m4mu->GetXaxis()->SetTitle("Invariant mass of 4muons (in GeV/c^2)"); h_m4_m4mu->GetYaxis()->SetTitle("Number of Events"); // These histograms are for 4 ELECTORON reconstruction with different combinations // First combination: 1234 // Mass Z12 TH1D *h_mZ12_4e = new TH1D("mZ12_4e", "mass of Z12", 75, 0., 150.); h_mZ12_4e->GetXaxis()->SetTitle("Invariant mass of dielectron (in GeV/c^2)"); h_mZ12_4e->GetYaxis()->SetTitle("Number of Events"); // Mass Z34 TH1D *h_mZ34_4e = new TH1D("mZ34_4e", "mass of Z34", 75, 0., 150.); h_mZ34_4e->GetXaxis()->SetTitle("Invariant mass of dielectron (in GeV/c^2)"); h_mZ34_4e->GetYaxis()->SetTitle("Number of Events"); // Second combination: 1324 // Mass Z13 TH1D *h_mZ13_4e = new TH1D("mZ13_4e", "mass of Z13", 75, 0., 150.); h_mZ13_4e->GetXaxis()->SetTitle("Invariant mass of dielectron (in GeV/c^2)"); h_mZ13_4e->GetYaxis()->SetTitle("Number of Events"); // Mass Z24 TH1D *h_mZ24_4e = new TH1D("mZ24_4e", "mass of Z24", 75, 0., 150.); h_mZ24_4e->GetXaxis()->SetTitle("Invariant mass of dielectron (in GeV/c^2)"); h_mZ24_4e->GetYaxis()->SetTitle("Number of Events"); // Third combination: 1423 // Mass Z14 TH1D *h_mZ14_4e = new TH1D("mZ14_4e", "mass of Z14", 75, 0., 150.); h_mZ14_4e->GetXaxis()->SetTitle("Invariant mass of dielectron (in GeV/c^2)"); h_mZ14_4e->GetYaxis()->SetTitle("Number of Events"); // Mass Z23 TH1D *h_mZ23_4e = new TH1D("mZ23_4e", "mass of Z23", 75, 0., 150.); h_mZ23_4e->GetXaxis()->SetTitle("Invariant mass of dielectron (in GeV/c^2)"); h_mZ23_4e->GetYaxis()->SetTitle("Number of Events"); // Mass Za: mass of Z closest to Z mass TH1D *h_mZa_4e = new TH1D("mZa_4e", "mass Za", 120, 0., 120.); h_mZa_4e->GetXaxis()->SetTitle("Invariant mass of dielectron (in GeV/c^2)"); h_mZa_4e->GetYaxis()->SetTitle("Number of Events"); // Mass Zb: mass of Z not closest to Z mass TH1D *h_mZb_4e = new TH1D("mZb_4e", "mass Zb", 120, 0., 120.); h_mZb_4e->GetXaxis()->SetTitle("Invariant mass of dielectron (in GeV/c^2)"); h_mZb_4e->GetYaxis()->SetTitle("Number of Events"); // 4electron mass (full mass range) TH1D *h_m4_m4e = new TH1D("mass4e_full", "mass of 4 electron", 300, 0., 900.); h_m4_m4e->GetXaxis()->SetTitle("Invariant mass of 4e (in GeV/c^2)"); h_m4_m4e->GetYaxis()->SetTitle("Number of Events"); // These histograms are for 2mu2e reconstruction with different combinations // Mass of Z to 2mu from 2mu2e TH1D *h_mZmu_2mu2e = new TH1D("massZmu_2mu2e", "mass Z2mu:2mu2e", 75, 0., 150.); h_mZmu_2mu2e->GetXaxis()->SetTitle("Invariant mass of Z1 (in GeV/c^2)"); h_mZmu_2mu2e->GetYaxis()->SetTitle("Number of Events"); // Mass of Z to 2e from 2mu2e TH1D *h_mZe_2mu2e = new TH1D("massZe_2mu2e", "mass Z2e:2mu2e", 75, 0., 150.); h_mZe_2mu2e->GetXaxis()->SetTitle("Invariant Mass of Z2 (in GeV/c^2)"); h_mZe_2mu2e->GetYaxis()->SetTitle("Number of Events"); // Mass Za: mass of Z1 closest to Z mass TH1D *h_mZa_2mu2e = new TH1D("mZa_2mu2e", "mass Z higher", 120, 0., 120.); h_mZa_2mu2e->GetXaxis()->SetTitle("Invariant mass of dimuon (in GeV/c^2)"); h_mZa_2mu2e->GetYaxis()->SetTitle("Number of Events"); // Mass Zb: mass of Z2 not closest to Z mass TH1D *h_mZb_2mu2e = new TH1D("mZb_2mu2e", "mass Z lower", 120, 0., 120.); h_mZb_2mu2e->GetXaxis()->SetTitle("Invariant mass of dimuon (in GeV/c^2)"); h_mZb_2mu2e->GetYaxis()->SetTitle("Number of Events"); // 2muons 2electrons mass spectrum full mass range TH1D *h_m4_m2mu2e = new TH1D("mass2mu2e_full", "mass of 2mu2e", 300, 0., 900.); h_m4_m2mu2e->GetXaxis()->SetTitle("Invariant mass of 2mu2e (in GeV/c^2)"); h_m4_m2mu2e->GetYaxis()->SetTitle("Number of Events"); // Initialize variables // select largest, init - pxZ12 = -9999.; pxZ34 = -9999.; pxZ13 = -9999.; pxZ24 = -9999.; pxZ14 = -9999.; pxZ23 = -9999.; pyZ12 = -9999.; pyZ34 = -9999.; pyZ13 = -9999.; pyZ24 = -9999.; pyZ14 = -9999.; pyZ23 = -9999.; pzZ12 = -9999.; pzZ34 = -9999.; pzZ13 = -9999.; pzZ24 = -9999.; pzZ14 = -9999.; pzZ23 = -9999.; mZ12 = -9999.; mZ34 = -9999.; mZ13 = -9999.; mZ24 = -9999.; mZ14 = -9999.; mZ23 = -9999.; dZ12 = 9999.; dZ34 = 9999.; dZ13 = 9999.; dZ24 = 9999.; dZ14 = 9999.; dZ23 = 9999.; // select smallest, init + dZc1 = 9999.; dZc2 = 9999.; dZc3 = 9999.; pxZa = -9999.; pyZa = -9999.; pzZa = -9999.; mZa = -9999.; pxZb = -9999.; pyZb = -9999.; pzZb = -9999.; mZb = -9999.; mass4mu = -9999.; mass4e = -9999.; mass2mu2e = -9999.; p4Za.SetPtEtaPhiM(0., 0., 0., 0.); p4Zb.SetPtEtaPhiM(0., 0., 0., 0.); p4H.SetPtEtaPhiM(0., 0., 0., 0.); // constant value squared muon mass, squared electron mass and Z mass mZ = 91.1876; // Loop over all events begins for(Int_t entry = 0; entry < numberOfEntries; ++entry) { // print event number in the file if(entry%10000==0) cout << "Reading Event " << entry << endl; // Load selected branches with data from specified event treeReader->ReadEntry(entry); electrons -> clear(); muons -> clear(); // Select leptons e||mu with pT > 5 || 7 GeV and |eta| < 2.5 || 2.4 CollectionFilter(*branchElectron, *electrons , 5.0 , 2.5); CollectionFilter(*branchMuon , *muons , 7.0 , 2.4); //sort the PT from highest to lowest // std::sort(ele.begin(), ele.end(), [] (const double &electron1, const double &electron2) {return electron1.PT < electron2.PT}); // sort Muon PT, High to low // std::sort(muons->begin(), muons->end(), [] (const double &mu1 , const double &mu2) { return mu1.PT < mu2.PT}); //=======================Zto2Muon Start=========================// if (muons->size() >= 2) { //define leading and subleading leptons muon1 = muons->at(0); muon2 = muons->at(1); //select event with opposite charge leptons if(muon1->Charge + muon2->Charge == 0) continue; //define muon pair invariant mass massMM = ((muon1->P4()) + (muon2->P4())).M(); h_mZ_2mu -> Fill(massMM); } //=========================Zto2e Start===========================// if(electrons->size() >= 2) continue; { //define leading and subleading leptons elec1 = electrons->at(0); elec2 = electrons->at(1); //select event with opposite charge leptons if(elec1->Charge + elec2->Charge == 0) continue; // Define electron pair invariant mass massEE = ((elec1->P4()) + (elec2->P4())).M(); h_mZ_2e -> Fill(massEE); } //==========================START ZZDto4Mu===============================// if (muons->size() >= 4) { muon1 = muons->at(0); muon2 = muons->at(1); muon3 = muons->at(2); muon4 = muons->at(3); // if (muon1->Charge + muon2->Charge + muon3->Charge + muon4->Charge == 0) if (muon1->Charge + muon2->Charge + muon3->Charge + muon4->Charge == 0) { //first combination: Combine muon 1234 if (muon1->Charge + muon2->Charge == 0) { pxZ12 = muon1->PT + muon2->PT; pyZ12 = muon1->Eta + muon2->Eta; pyZ12 = muon1->Phi + muon2->Phi; /*eZ12 = (sqrt(muon1.P() * muon1.P() + sqm1)) + (sqrt(muon2.P() * muon2.P() + sqm1)); pxZ12 = muon1.Px() + muon2.Px(); pyZ12 = muon1.Py() + muon2.Py(); pzZ12 = muon1.Pz() + muon2.Pz(); */ if (muon3->Charge + muon4->Charge == 0) { pxZ34 = muon3->PT + muon4->PT; pyZ34 = muon3->Eta + muon4->Eta; pyZ34 = muon3->Phi + muon4->Phi; /*eZ34 = (sqrt(muon3.P() * muon3.P() + sqm1)) + (sqrt(muon4.P() * muon4.P() + sqm1)); pxZ34 = muon3.Px() + muon4.Px(); pyZ34 = muon3.Py() + muon4.Py(); pzZ34 = muon3.Pz() + muon4.Pz(); // Calculate p4 pZ12 = sqrt((pxZ12 * pxZ12) + (pyZ12 * pyZ12) + (pzZ12 * pzZ12)); pZ34 = sqrt((pxZ34 * pxZ34) + (pyZ34 * pyZ34) + (pzZ34 * pzZ34)); pTZ12 = sqrt((pxZ12 * pxZ12) + (pyZ12 * pyZ12)); pTZ34 = sqrt((pxZ34 * pxZ34) + (pyZ34 * pyZ34)); mZ12 = sqrt((eZ12 * eZ12) - (pZ12 * pZ12)); mZ34 = sqrt((eZ34 * eZ34) - (pZ34 * pZ34)); */ mZ12 = ((muon1->P4()) + (muon2->P4())).M(); mZ34 = ((muon3->P4()) + (muon4->P4())).M(); if (mZ12 > 0.) h_mZ12_4mu->Fill(mZ12); if (mZ34 > 0.) h_mZ34_4mu->Fill(mZ34); } } //close first combination mu1234 dZ12 = std::abs( mZ12 - mZ ); dZ34 = std::abs( mZ34 - mZ ); // take the smallest difference between mass // to use for 4muon combination dZc1 = (dZ12 < dZ34) ? dZ12 : dZ34; // Second combination: Combine muon 1324 if (muon1->Charge + muon3->Charge == 0) { pxZ13 = muon1->PT + muon3->PT; pyZ13 = muon1->Eta + muon3->Eta; pyZ13 = muon1->Phi + muon3->Phi; /*eZ13 = (sqrt(muon1.P() * muon1.P() + sqm1)) + (sqrt(muon3.P() * muon3.P() + sqm1)); pxZ13 = muon1.Px() + muon3.Px(); pyZ13 = muon1.Py() + muon3.Py(); pzZ13 = muon1.Pz() + muon3.Pz();*/ if (muon2->Charge + muon4->Charge == 0) { pxZ24 = muon2->PT + muon4->PT; pyZ24 = muon2->Eta + muon4->Eta; pzZ24 = muon2->Phi + muon2->Phi; /*eZ24 = (sqrt(muon2.P() * muon2.P() + sqm1)) + (sqrt(muon4.P() * muon4.P() + sqm1)); pxZ24 = muon2.Px() + muon4.Px(); pyZ24 = muon2.Py() + muon4.Py(); pzZ24 = muon2.Pz() + muon4.Pz(); // Calculate p4 pZ13 = sqrt((pxZ13 * pxZ13) + (pyZ13 * pyZ13) + (pzZ13 * pzZ13)); pZ24 = sqrt((pxZ24 * pxZ24) + (pyZ24 * pyZ24) + (pzZ24 * pzZ24)); pTZ13 = sqrt((pxZ13 * pxZ13) + (pyZ13 * pyZ13)); pTZ24 = sqrt((pxZ24 * pxZ24) + (pyZ24 * pyZ24)); mZ13 = sqrt((eZ13 * eZ13) - (pZ13 * pZ13)); mZ24 = sqrt((eZ24 * eZ24) - (pZ24 * pZ24)); */ mZ13 = ((muon1->P4()) + (muon3->P4())).M(); mZ24 = ((muon2->P4()) + (muon4->P4())).M(); if (mZ13 > 0.) h_mZ13_4mu->Fill(mZ13); if (mZ24 > 0.) h_mZ24_4mu->Fill(mZ24); } } //close second combination mu1324 dZ13 = std::abs( mZ13 - mZ ); dZ24 = std::abs( mZ24 - mZ ); dZc2 = (dZ13 < dZ24) ? dZ13 : dZ24; // Third combination: Combine muon 1423 if (muon1->Charge + muon4->Charge == 0) { pxZ14 = muon1->PT + muon4->PT; pyZ14 = muon1->Eta + muon4->Eta; pzZ14 = muon1->Phi + muon4->Phi; /*eZ14 = (sqrt(muon1.P() * muon1.P() + sqm1)) + (sqrt(muon4.P() * muon4.P() + sqm1)); pxZ14 = pxZ14 = muon1.Px() + muon4.Px(); pyZ14 = muon1.Py() + muon4.Py(); pzZ14 = muon1.Pz() + muon4.Pz();*/ if (muon2->Charge + muon3->Charge == 0) { pxZ23 = muon2->PT + muon3->PT; pyZ23 = muon2->Eta + muon3->Eta; pzZ23 = muon2->Phi + muon3->Phi; /*eZ23 = sqrt((muon2.P() * muon2.P() + sqm1)) + (sqrt(muon3.P() * muon3.P() + sqm1)); pxZ23 = muon2.Px() + muon3.Px(); pyZ23 = muon2.Py() + muon3.Py(); pzZ23 = muon2.Pz() + muon3.Pz(); // Calculate p4 pZ14 = sqrt((pxZ14 * pxZ14) + (pyZ14 * pyZ14) + (pzZ14 * pzZ14)); pZ23 = sqrt((pxZ23 * pxZ23) + (pyZ23 * pyZ23) + (pzZ23 * pzZ23)); pTZ14 = sqrt((pxZ14 * pxZ14) + (pyZ14 * pyZ14)); pTZ23 = sqrt((pxZ23 * pxZ23) + (pyZ23 * pyZ23)); mZ14 = sqrt((eZ14 * eZ14) - (pZ14 * pZ14)); mZ23 = sqrt((eZ23 * eZ23) - (pZ23 * pZ23)); */ mZ14 = ((muon1->P4()) + (muon4->P4())).M(); mZ23 = ((muon2->P4()) + (muon3->P4())).M(); if (mZ14 > 0.) h_mZ14_4mu->Fill(mZ14); if (mZ23 > 0.) h_mZ23_4mu->Fill(mZ23); } } // close thrid combination mu1423 dZ14 = std::abs( mZ14 - mZ ); dZ23 = std::abs( mZ23 - mZ ); dZc3 = (dZ14 < dZ23) ? dZ14 : dZ23; bool ptZadaug = false; if (dZc1 < dZc2 && dZc1 < dZc3) { if (dZ12 < dZ34) { //eZa = eZ12; pxZa = pxZ12; pyZa = pyZ12; pzZa = pzZ12; // pTZa = pTZ12; mZa = mZ12; if (muon1->PT > 20. and muon2->PT > 10.) ptZadaug = true; //eZb = eZ34; pxZb = pxZ34; pyZb = pyZ34; pzZb = pzZ34; //pTZb = pTZ34; mZb = mZ34; } else { //eZa = eZ34; pxZa = pxZ34; pyZa = pyZ34; pzZa = pzZ34; //pTZa = pTZ34; mZa = mZ34; if (muon3->PT > 20. and muon4->PT > 10.) ptZadaug = true; //eZb = eZ12; pxZb = pxZ12; pyZb = pyZ12; pzZb = pzZ12; //pTZb = pTZ12; mZb = mZ12; } } // close (dZc1 < dZc2 && dZc1 < dZc3) else if (dZc2 < dZc1 && dZc2 < dZc3) { if (dZ13 < dZ24) { //eZa = eZ13; pxZa = pxZ13; pyZa = pyZ13; pzZa = pzZ13; //pTZa = pTZ13; mZa = mZ13; if (muon1->PT > 20. and muon3->PT > 10.) ptZadaug = true; //eZb = eZ24; pxZb = pxZ24; pyZb = pyZ24; pzZb = pzZ24; // pTZb = pTZ24; mZb = mZ24; } else { //eZa = eZ24; pxZa = pxZ24; pyZa = pyZ24; pzZa = pzZ24; //pTZa = pTZ24; mZa = mZ24; if (muon2->PT > 20. and muon4->PT > 10.) ptZadaug = true; //eZb = eZ13; pxZb = pxZ13; pyZb = pyZ13; pzZb = pzZ13; //pTZb = pTZ13; mZb = mZ13; } } // close (dZc2 < dZc1 && dZc2 < dZc3) else if (dZc3 < dZc1 && dZc3 < dZc2) { if (dZ14 < dZ23) { //eZa = eZ14; pxZa = pxZ14; pyZa = pyZ14; pzZa = pzZ14; // pTZa = pTZ14; mZa = mZ14; if (muon1->PT > 20. and muon4->PT > 10.) ptZadaug = true; //eZb = eZ23; pxZb = pxZ23; pyZb = pyZ23; pzZb = pzZ23; // pTZb = pTZ23; mZb = mZ23; } else { //eZa = eZ23; pxZa = pxZ23; pyZa = pyZ23; pzZa = pzZ23; // pTZa = pTZ23; mZa = mZ23; if (muon2->PT > 20. and muon3->PT > 10.) ptZadaug = true; //eZb = eZ14; pxZb = pxZ14; pyZb = pyZ14; pzZb = pzZ14; //pTZb = pTZ14; mZb = mZ14; } } // close (dZc3 < dZc1 && dZc3 < dZc2) if (ptZadaug) { if (mZa > 40. && mZa < 120.) { if (mZb > 12. && mZb < 120.){ h_mZa_4mu->Fill(mZa); h_mZb_4mu->Fill(mZb); // 4 vector p4Za.SetPtEtaPhiM(pxZa, pyZa, pzZa, mZa); p4Zb.SetPtEtaPhiM(pxZb, pyZb, pzZb, mZb); p4H = p4Za + p4Zb; mass4mu = p4H.M(); /* pt_4mu = p4H.Pt(); eta_4mu = p4H.Eta(); phi_4mu = p4H.Phi(); //px4mu = p4H.Px(); //py4mu = p4H.Py(); // pz4mu = p4H.Pz(); // E4mu = p4H.E(); pt_mu1 = muon1->PT; pt_mu2 = muon2->PT; pt_mu3 = muon3->PT; pt_mu4 = muon4->PT; eta_mu1 = muon1->Eta; eta_mu2 = muon2->Eta; eta_mu3 = muon3->Eta; eta_mu4 = muon4->Eta; phi_mu1 = muon1->Phi; phi_mu2 = muon2->Phi; phi_mu3 = muon3->Phi; phi_mu4 = muon4->Phi; cas_mu1 = muon1->Charge; cas_mu2 = muon2->Charge; cas_mu3 = muon3->Charge; cas_mu4 = muon4->Charge; //px_mu1 = muon1.Px(); // px_mu2 = muon2.Px(); // px_mu3 = muon3.Px(); // px_mu4 = muon4.Px(); //py_mu1 = muon1.Py(); // py_mu2 = muon2.Py(); // py_mu3 = muon3.Py(); // py_mu4 = muon4.Py(); //pz_mu1 = muon1.Pz(); //pz_mu2 = muon2.Pz(); //pz_mu3 = muon3.Pz(); //pz_mu4 = muon4.Pz(); // E_mu1 = sqrt((muon1.P() * muon1.P()) + sqm1); // E_mu2 = sqrt((muon2.P() * muon2.P()) + sqm1); /// E_mu3 = sqrt((muon3.P() * muon3.P()) + sqm1); // E_mu4 = sqrt((muon4.P() * muon4.P()) + sqm1); */ if (mass4mu > 70.) { h_m4_m4mu->Fill(mass4mu); } } // close zb } //close za } //end ptZadaug }//end 4mu charge } // close All 4mu //****************************************END 4MU************************************// //===================================START ZZDTo4e==================================// if (electrons->size() >= 4 ) { elec1 = electrons->at(0); elec2 = electrons->at(1); elec3 = electrons->at(2); elec4 = electrons->at(3); if (elec1->Charge + elec2->Charge + elec3->Charge + elec4->Charge == 0) { //FIRST combination : combine elec 1234 if (elec1->Charge + elec2->Charge == 0) { pxZ12 = elec1->PT + elec2->PT; pyZ12 = elec1->Eta + elec2->Eta; pzZ12 = elec1->Phi + elec2->Phi; /*eZ12 = (sqrt(elec1.P() * elec1.P() + sqme)) + (sqrt(elec2.P() * elec2.P() + sqme)); pxZ12 = elec1.Px() + elec2.Px(); pyZ12 = elec1.Py() + elec2.Py(); pzZ12 = elec1.Pz() + elec2.Pz();*/ if (elec3->Charge + elec4->Charge == 0) { pxZ34 = elec3->PT + elec4->PT; pyZ34 = elec3->Eta + elec4->Eta; pzZ34 = elec3->Phi + elec4->Phi; /*eZ34 = (sqrt(elec3.P() * elec3.P() + sqme)) + (sqrt(elec4.P() * elec4.P() + sqme)); pxZ34 = elec3.Px() + elec4.Px(); pyZ34 = elec3.Py() + elec4.Py(); pzZ34 = elec3.Pz() + elec4.Pz(); // Calculate the momentum and invariant mass of elec 1 and 2, elec 3 and 4 pZ12 = sqrt((pxZ12 * pxZ12) + (pyZ12 * pyZ12) + (pzZ12 * pzZ12)); pZ34 = sqrt((pxZ34 * pxZ34) + (pyZ34 * pyZ34) + (pzZ34 * pzZ34)); pTZ12 = sqrt((pxZ12 * pxZ12) + (pyZ12 * pyZ12)); pTZ34 = sqrt((pxZ34 * pxZ34) + (pyZ34 * pyZ34)); mZ12 = sqrt((eZ12 * eZ12) - (pZ12 * pZ12)); mZ34 = sqrt((eZ34 * eZ34) - (pZ34 * pZ34));*/ if (mZ12 > 0.) h_mZ12_4e->Fill(mZ12); if (mZ34 > 0.) h_mZ34_4e->Fill(mZ34); } } // close combination 1, elec 1234 dZ12 = std::abs( mZ12 - mZ ); dZ34 = std::abs( mZ34 - mZ ); // take the smallest diff between mass to use for 4electron combination dZc1 = (dZ12 < dZ34) ? dZ12 : dZ34; // Second combination: Combine elec 1324 if (elec1->Charge + elec3->Charge == 0) { pxZ13 = elec1->PT + elec3->PT; pyZ13 = elec1->Eta + elec3->Eta; pzZ13 = elec1->Phi + elec3->Phi; /*eZ13 = (sqrt(elec1.P() * elec1.P() + sqme)) + (sqrt(elec3.P() * elec3.P() + sqme)); pxZ13 = elec1.Px() + elec3.Px(); pyZ13 = elec1.Py() + elec3.Py(); pzZ13 = elec1.Pz() + elec3.Pz();*/ if (elec2->Charge + elec4->Charge == 0) { pxZ24 = elec2->PT + elec4->PT; pyZ24 = elec2->Eta + elec4->Eta; pzZ24 = elec2->Phi + elec4->Phi; /*eZ24 = (sqrt(elec2.P() * elec2.P() + sqme)) + (sqrt(elec4.P() * elec4.P() + sqme)); pxZ24 = elec2.Px() + elec4.Px(); pyZ24 = elec2.Py() + elec4.Py(); pzZ24 = elec2.Pz() + elec4.Pz(); // Calculate the momentum and invariant mass of elec 1 and 3, elec 2 and 4 pZ13 = sqrt((pxZ13 * pxZ13) + (pyZ13 * pyZ13) + (pzZ13 * pzZ13)); pZ24 = sqrt((pxZ24 * pxZ24) + (pyZ24 * pyZ24) + (pzZ24 * pzZ24)); pTZ13 = sqrt((pxZ13 * pxZ13) + (pyZ13 * pyZ13) + (pzZ13 * pzZ13)); pTZ24 = sqrt((pxZ24 * pxZ24) + (pyZ24 * pyZ24) + (pzZ24 * pzZ24)); mZ13 = sqrt((eZ13 * eZ13) - (pZ13 * pZ13)); mZ24 = sqrt((eZ24 * eZ24) - (pZ24 * pZ24));*/ if (mZ13 > 0.) h_mZ13_4e->Fill(mZ13); if (mZ24 > 0.) h_mZ24_4e->Fill(mZ24); } } //close second combination 1324 dZ13 = std::abs( mZ13 - mZ ); dZ24 = std::abs( mZ24 - mZ ); dZc2 = (dZ13 < dZ24) ? dZ13 : dZ24; // Third combination: Combine elec 1423 if (elec1->Charge + elec4->Charge == 0) { pxZ14 = elec1->PT + elec4->PT; pyZ14 = elec1->Eta + elec4->Eta; pzZ14 = elec1->Phi + elec4->Phi; /*eZ14 = (sqrt(elec1.P() * elec1.P() + sqme)) + (sqrt(elec4.P() * elec4.P() + sqme)); pxZ14 = elec1.Px() + elec4.Px(); pyZ14 = elec1.Py() + elec4.Py(); pzZ14 = elec1.Pz() + elec4.Pz();*/ if (elec2->Charge + elec3->Charge == 0) { pxZ23 = elec2->PT + elec3->PT; pyZ23 = elec2->Eta + elec3->Eta; pzZ23 = elec1->Phi + elec3->Phi; /*eZ23 = sqrt((elec2.P() * elec2.P() + sqme)) + (sqrt(elec3.P() * elec3.P() + sqme)); pxZ23 = elec2.Px() + elec3.Px(); pyZ23 = elec2.Py() + elec3.Py(); pzZ23 = elec2.Pz() + elec3.Pz(); // Calculate the momentum and invariant mass of elec 1 and 4, elec 2 and 3 pZ14 = sqrt((pxZ14 * pxZ14) + (pyZ14 * pyZ14) + (pzZ14 * pzZ14)); pZ23 = sqrt((pxZ23 * pxZ23) + (pyZ23 * pyZ23) + (pzZ23 * pzZ23)); pTZ14 = sqrt((pxZ14 * pxZ14) + (pyZ14 * pyZ14) + (pzZ14 * pzZ14)); pTZ23 = sqrt((pxZ23 * pxZ23) + (pyZ23 * pyZ23) + (pzZ23 * pzZ23)); mZ14 = sqrt((eZ14 * eZ14) - (pZ14 * pZ14)); mZ23 = sqrt((eZ23 * eZ23) - (pZ23 * pZ23));*/ if (mZ14 > 0.) h_mZ14_4e->Fill(mZ14); if (mZ23 > 0.) h_mZ23_4e->Fill(mZ23); } } // close third combination 1423 dZ14 = std::abs( mZ14 - mZ ); dZ23 = std::abs( mZ23 - mZ ); dZc3 = (dZ14 < dZ23) ? dZ14 : dZ23; bool ptZadaug = false; // Now whichever have the smallest diff is considered the best comb. if (dZc1 < dZc2 && dZc1 < dZc3) { if (dZ12 < dZ34) { //eZa = eZ12; pxZa = pxZ12; pyZa = pyZ12; pzZa = pzZ12; //pTZa = pTZ12; mZa = mZ12; if (elec1->PT > 20. and elec2->PT > 10.) ptZadaug = true; // eZb = eZ34; pxZb = pxZ34; pyZb = pyZ34; pzZb = pzZ34; //pTZb = pTZ34; mZb = mZ34; } else { //eZa = eZ34; pxZa = pxZ34; pyZa = pyZ34; pzZa = pzZ34; // pTZa = pTZ34; mZa = mZ34; if (elec3->PT > 20. and elec4->PT > 10.) ptZadaug = true; //eZb = eZ12; pxZb = pxZ12; pyZb = pyZ12; pzZb = pzZ12; // pTZb = pTZ12; mZb = mZ12; } } // close (dZc1 < dZc2 && dZc1 < dZc3) else if (dZc2 < dZc1 && dZc2 < dZc3) { if (dZ13 < dZ24) { //eZa = eZ13; pxZa = pxZ13; pyZa = pyZ13; pzZa = pzZ13; // pTZa = pTZ13; mZa = mZ13; if (elec1->PT > 20. and elec3->PT > 10.) ptZadaug = true; //eZb = eZ24; pxZb = pxZ24; pyZb = pyZ24; pzZb = pzZ24; // pTZb = pTZ24; mZb = mZ24; } else { // eZa = eZ24; pxZa = pxZ24; pyZa = pyZ24; pzZa = pzZ24; // pTZa = pTZ24; mZa = mZ24; if (elec2->PT > 20. and elec4->PT > 10.) ptZadaug = true; // eZb = eZ13; pxZb = pxZ13; pyZb = pyZ13; pzZb = pzZ13; // pTZb = pTZ13; mZb = mZ13; } } //close (dZc2 < dZc1 && dZc2 < dZc3) else if (dZc3 < dZc1 && dZc3 < dZc2) { if (dZ14 < dZ23) { // eZa = eZ14; pxZa = pxZ14; pyZa = pyZ14; pzZa = pzZ14; // pTZa = pTZ14; mZa = mZ14; if (elec1->PT > 20. and elec4->PT > 10.) ptZadaug = true; //eZb = eZ23; pxZb = pxZ23; pyZb = pyZ23; pzZb = pzZ23; // pTZb = pTZ23; mZb = mZ23; } else { // eZa = eZ23; pxZa = pxZ23; pyZa = pyZ23; pzZa = pzZ23; //pTZa = pTZ23; mZa = mZ23; if (elec2->PT > 20. and elec3->PT > 10.) ptZadaug = true; eZb = eZ14; pxZb = pxZ14; pyZb = pyZ14; pzZb = pzZ14; pTZb = pTZ14; mZb = mZ14; } } //close (dZc3 < dZc1 && dZc3 < dZc2) if (ptZadaug) { if (mZa > 40. && mZa < 120.) { if (mZb > 12. && mZb < 120.) { h_mZa_4e->Fill(mZa); h_mZb_4e->Fill(mZb); // Calculate 4 elec //p4Za.SetPxPyPzE(pxZa, pyZa, pzZa, eZa); //p4Zb.SetPxPyPzE(pxZb, pyZb, pzZb, eZb); p4Za.SetPtEtaPhiM(pxZa, pyZa, pzZa, mZa); p4Zb.SetPtEtaPhiM(pxZb, pyZb, pzZb, mZb); p4H = p4Za + p4Zb; mass4e = p4H.M(); /* pt_4e = p4H.PT(); eta_4e = p4H.Eta(); phi_4e = p4H.Phi(); px4e = p4H.Px(); py4e = p4H.Py(); pz4e = p4H.Pz(); E4e = p4H.E(); pt_e1 = elec1.PT(); pt_e2 = elec2.PT(); pt_e3 = elec3.PT(); pt_e4 = elec4.PT(); eta_e1 = elec1.Eta(); eta_e2 = elec2.Eta(); eta_e3 = elec3.Eta(); eta_e4 = elec4.Eta(); phi_e1 = elec1.Phi(); phi_e2 = elec2.Phi(); phi_e3 = elec3.Phi(); phi_e4 = elec4.Phi(); cas_e1 = elec1.Charge(); cas_e2 = elec2.Charge(); cas_e3 = elec3.Charge(); cas_e4 = elec4.Charge(); px_e1 = elec1.Px(); px_e2 = elec2.Px(); px_e3 = elec3.Px(); px_e4 = elec4.Px(); py_e1 = elec1.Py(); py_e2 = elec2.Py(); py_e3 = elec3.Py(); py_e4 = elec4.Py(); pz_e1 = elec1.Pz(); pz_e2 = elec2.Pz(); pz_e3 = elec3.Pz(); pz_e4 = elec4.Pz(); E_e1 = elec1.E(); E_e2 = elec2.E(); E_e3 = elec3.E(); E_e4 = elec4.E(); */ if (mass4e > 70.) { h_m4_m4e->Fill(mass4e); } } //end zb } //end za } // end ptZadaug }//close 4e charge } //close 4e //************************************** END ZZDTo4e ************************************// //=================================== START ZZD To2mu2e =================================// if (electrons->size() >= 2 && muons->size() >= 2 ) { elec1 = electrons->at(0); elec2 = electrons->at(1); muon1 = muons->at(0); muon2 = muons->at(1); if (muon1->Charge + muon2->Charge + elec1->Charge + elec2->Charge == 0) { // For case 2mu2e, there is only 1 combination if (muon1->Charge + muon2->Charge == 0) { pxZ12 = muon1->PT + muon2->PT; pyZ12 = muon1->Eta + muon2->Eta; pzZ12 = muon1->Phi + muon2->Phi; /*eZ12 = (sqrt(muon1.P() * muon1.P() + sqm1)) + (sqrt(muon2.P() * muon2.P() + sqm1)); pxZ12 = muon1.Px() + muon2.Px(); pyZ12 = muon1.Py() + muon2.Py(); pzZ12 = muon1.Pz() + muon2.Pz();*/ if (elec1->Charge + elec2->Charge == 0) { pxZ12 = elec1->PT + elec2->PT; pyZ12 = elec1->Eta + elec2->Eta; pzZ12 = elec1->Phi + elec2->Phi; /*eZ34 = (sqrt(elec1.P() * elec1.P() + sqme)) + (sqrt(elec2.P() * elec2.P() + sqme)); pxZ34 = elec1.Px() + elec2.Px(); pyZ34 = elec1.Py() + elec2.Py(); pzZ34 = elec1.Pz() + elec2.Pz(); // Calculate the momentum and invariant mass of muon 12, elec 34 pZ12 = sqrt((pxZ12 * pxZ12) + (pyZ12 * pyZ12) + (pzZ12 * pzZ12)); pZ34 = sqrt((pxZ34 * pxZ34) + (pyZ34 * pyZ34) + (pzZ34 * pzZ34)); pTZ12 = sqrt((pxZ12 * pxZ12) + (pyZ12 * pyZ12)); pTZ34 = sqrt((pxZ34 * pxZ34) + (pyZ34 * pyZ34)); mZ12 = sqrt((eZ12 * eZ12) - (pZ12 * pZ12)); mZ34 = sqrt((eZ34 * eZ34) - (pZ34 * pZ34));*/ if (mZ12 > 0.) h_mZmu_2mu2e->Fill(mZ12); if (mZ34 > 0.) h_mZe_2mu2e->Fill(mZ34); } } dZ12 = std::abs(mZ12 - mZ); // mu dZ34 = std::abs(mZ34 - mZ); // e bool ptZadaug = false; if (dZ12 < dZ34) { //eZa = eZ12; pxZa = pxZ12; pyZa = pyZ12; pzZa = pzZ12; // pTZa = pTZ12; mZa = mZ12; if (muon1->PT > 20. and muon2->PT > 10.) ptZadaug = true; //eZb = eZ34; pxZb = pxZ34; pyZb = pyZ34; pzZb = pzZ34; //pTZb = pTZ34; mZb = mZ34; } else { //eZa = eZ34; pxZa = pxZ34; pyZa = pyZ34; pzZa = pzZ34; //pTZa = pTZ34; mZa = mZ34; if (elec1->PT > 20. and elec2->PT > 10.) ptZadaug = true; // eZb = eZ12; pxZb = pxZ12; pyZb = pyZ12; pzZb = pzZ12; // pTZb = pTZ12; mZb = mZ12; } if (ptZadaug) { if (mZa > 40. && mZa < 120.) { if (mZb > 12. && mZb < 120.) { h_mZa_2mu2e->Fill(mZa); h_mZb_2mu2e->Fill(mZb); // Now combine these 2 muons and 2 electrons // Calculate 4 lepton: 2muons 2electrons p4Za.SetPtEtaPhiM(pxZa, pyZa, pzZa, mZa); p4Zb.SetPtEtaPhiM(pxZb, pyZb, pzZb, mZb); p4H = p4Za + p4Zb; mass2mu2e = p4H.M(); /* pt_2mu2e = p4H.PT(); eta_2mu2e = p4H.Eta(); phi_2mu2e = p4H.Phi(); px2mu2e = p4H.Px(); py2mu2e = p4H.Py(); pz2mu2e = p4H.Pz(); E2mu2e = p4H.E(); pt_2mu1 = muon1.PT(); pt_2mu2 = muon2.PT(); pt_2e1 = elec1.PT(); pt_2e2 = elec2.PT(); eta_2mu1 = muon1.Eta(); eta_2mu2 = muon2.Eta(); eta_2e1 = elec1.Eta(); eta_2e2 = elec2.Eta(); phi_2mu1 = muon1.Phi(); phi_2mu2 = muon2.Phi(); phi_2e1 = elec1.Phi(); phi_2e2 = elec2.Phi(); cas_2mu1 = muon1.Charge(); cas_2mu2 = muon2.Charge(); cas_2e1 = elec1.Charge(); cas_2e2 = elec2.Charge(); px_2mu1 = muon1.Px(); px_2mu2 = muon2.Px(); px_2e1 = elec1.Px(); px_2e2 = elec2.Px(); py_2mu1 = muon1.Py(); py_2mu2 = muon2.Py(); py_2e1 = elec1.Py(); py_2e2 = elec2.Py(); pz_2mu1 = muon1.Pz(); pz_2mu2 = muon2.Pz(); pz_2e1 = elec1.Pz(); pz_2e2 = elec2.Pz(); E_2mu1 = sqrt((muon1.P() * muon1.P()) + sqm1); E_2mu2 = sqrt((muon2.P() * muon2.P()) + sqm1); E_2e1 = elec1.E(); E_2e2 = elec2.E();*/ if (mass2mu2e > 70.) { //h_m1_m2mu2e->Fill(mass2mu2e); // h_m2_m2mu2e->Fill(mass2mu2e); // h_m3_m2mu2e->Fill(mass2mu2e); h_m4_m2mu2e->Fill(mass2mu2e); } } //close zb }//close za }//close ptZadaug }//close charge } // close 2mu2e //*************************************ZZdTo2e2mu end********************************// } // END EVENT LOOP //CREATE CANVAS TCanvas *c = new TCanvas("C1","4LANALYZER",200,200,600,600); // CREATE AND OPEN NEW ROOT file to store info TFile fout("Audiya4lAnalyzer.root", "RECREATE"); //Write the histograms //dimuon pair h_mZ_2mu->Write(); h_mZ_2e->Write(); //4mu h_mZ12_4mu->Write(); h_mZ34_4mu->Write(); h_mZ13_4mu->Write(); h_mZ24_4mu->Write(); h_mZ14_4mu->Write(); h_mZ23_4mu->Write(); h_mZa_4mu->Write(); h_mZb_4mu->Write(); h_m4_m4mu->Write(); //4e h_mZ12_4e->Write(); h_mZ34_4e->Write(); h_mZ13_4e->Write(); h_mZ24_4e->Write(); h_mZ14_4e->Write(); h_mZ23_4e->Write(); h_mZa_4e->Write(); h_mZb_4e->Write(); h_m4_m4e->Write(); //2e2mu h_mZmu_2mu2e->Write(); h_mZe_2mu2e->Write(); h_mZa_2mu2e->Write(); h_mZb_2mu2e->Write(); h_m4_m2mu2e->Write(); // Clost the ROOT File fout.Close(); //Draw the histogram //dimuon pair 2 h_mZ_2mu->Draw(); h_mZ_2e->Draw(); //4mu 9 h_mZ12_4mu->Draw(); h_mZ34_4mu->Draw(); h_mZ13_4mu->Draw(); h_mZ24_4mu->Draw(); h_mZ14_4mu->Draw(); h_mZ23_4mu->Draw(); h_mZa_4mu->Draw(); h_mZb_4mu->Draw(); h_m4_m4mu->Draw(); //4e 9 h_mZ12_4e->Draw(); h_mZ34_4e->Draw(); h_mZ13_4e->Draw(); h_mZ24_4e->Draw(); h_mZ14_4e->Draw(); h_mZ23_4e->Draw(); h_mZa_4e->Draw(); h_mZb_4e->Draw(); h_m4_m4e->Draw(); //2e2mu 5 h_mZmu_2mu2e->Draw(); h_mZe_2mu2e->Draw(); h_mZa_2mu2e->Draw(); h_mZb_2mu2e->Draw(); h_m4_m2mu2e->Draw(); } // END VOID