当前位置: 首页>>代码示例>>C++>>正文


C++ TLorentzVector::Rapidity方法代码示例

本文整理汇总了C++中TLorentzVector::Rapidity方法的典型用法代码示例。如果您正苦于以下问题:C++ TLorentzVector::Rapidity方法的具体用法?C++ TLorentzVector::Rapidity怎么用?C++ TLorentzVector::Rapidity使用的例子?那么, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在TLorentzVector的用法示例。


在下文中一共展示了TLorentzVector::Rapidity方法的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。

示例1: AnalyseParticles

void AnalyseParticles(LHEF::Reader *reader, ExRootTreeBranch *branch)
{
  const LHEF::HEPEUP &hepeup = reader->hepeup;

  Int_t particle;
  Double_t signPz, cosTheta;

  TLorentzVector momentum;

  TRootLHEFParticle *element;

  for(particle = 0; particle < hepeup.NUP; ++particle)
  {
    element = (TRootLHEFParticle*) branch->NewEntry();

    element->PID = hepeup.IDUP[particle];
    element->Status = hepeup.ISTUP[particle];
    element->Mother1 = hepeup.MOTHUP[particle].first;
    element->Mother2 = hepeup.MOTHUP[particle].second;
    element->ColorLine1 = hepeup.ICOLUP[particle].first;
    element->ColorLine2 = hepeup.ICOLUP[particle].second;
    element->Px = hepeup.PUP[particle][0];
    element->Py = hepeup.PUP[particle][1];
    element->Pz = hepeup.PUP[particle][2];
    element->E = hepeup.PUP[particle][3];
    element->M = hepeup.PUP[particle][4];

    momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);

    cosTheta = TMath::Abs(momentum.CosTheta());
    signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;

    element->PT = momentum.Perp();
    element->Eta = (cosTheta == 1.0 ? signPz*999.9 : momentum.Eta());
    element->Phi = (cosTheta == 1.0 ? signPz*999.9 : momentum.Rapidity());
    element->Rapidity = (cosTheta == 1.0 ? signPz*999.9 : momentum.Rapidity());

    element->LifeTime = hepeup.VTIMUP[particle];
    element->Spin = hepeup.SPINUP[particle];
  }
}
开发者ID:brovercleveland,项目名称:HZG_Analyzer,代码行数:41,代码来源:ExRootLHEFConverter.cpp

示例2: loop


//.........这里部分代码省略.........
			}
		      Btypesize[tidx]++;
		    }
		}
	      if(t!=3)
		{
		  if(ptflag>=0) Bmaxpt[ptflag] = true;
		  if(probflag>=0) Bmaxprob[probflag] = true;
		  if(tktkflag>=0) Bbesttktkmass[tktkflag] = true;
		  if(ptMatchedflag>=0) BmaxptMatched[ptMatchedflag] = true;
		  if(probMatchedflag>=0) BmaxprobMatched[probMatchedflag] = true;
		  if(tktkMatchedflag>=0) BbesttktkmassMatched[tktkMatchedflag] = true;
		}
	      if(t==0)      nt0->Fill();
	      else if(t==1) nt1->Fill();
	      else if(t==2) nt2->Fill();
	      else if(t==4) nt3->Fill();
	      else if(t==5) nt5->Fill();
	      else if(t==6) nt6->Fill();
	    }
	}

      ntHlt->Fill();
      if(isPbPb) ntHi->Fill();

      if(!REAL)
	{
	  Int_t gt=0,sigtype=0;
	  Int_t gsize=0;
	  Gsize = 0;
	  for(int j=0;j<GenInfo_size;j++)
	    {
	      if((TMath::Abs(GenInfo_pdgId[j])!=BPLUS_PDGID&&TMath::Abs(GenInfo_pdgId[j])!=BZERO_PDGID&&TMath::Abs(GenInfo_pdgId[j])!=BSUBS_PDGID) && gskim) continue;
	      Gsize = gsize+1;
	      Gpt[gsize] = GenInfo_pt[j];
	      Geta[gsize] = GenInfo_eta[j];
	      Gphi[gsize] = GenInfo_phi[j];
	      GpdgId[gsize] = GenInfo_pdgId[j];
	      bGen->SetPtEtaPhiM(GenInfo_pt[j],GenInfo_eta[j],GenInfo_phi[j],GenInfo_mass[j]);
	      Gy[gsize] = bGen->Rapidity();
	      sigtype=0;
	      for(gt=1;gt<8;gt++)
		{
		  if(signalGen(gt,j))
		    {
		      sigtype=gt;
		      break;
		    }
		}
	      GisSignal[gsize] = sigtype;
	      Gmu1pt[gsize] = -1;
	      Gmu1eta[gsize] = -20;
	      Gmu1phi[gsize] = -20;
	      Gmu1p[gsize] = -1;
	      Gmu2pt[gsize] = -1;
	      Gmu2eta[gsize] = -20;
	      Gmu2phi[gsize] = -20;
	      Gmu2p[gsize] = -1;
	      Gtk1pt[gsize] = -1;
	      Gtk1eta[gsize] = -20;
	      Gtk1phi[gsize] = -20;
	      Gtk2pt[gsize] = -1;
	      Gtk2eta[gsize] = -20;
	      Gtk2phi[gsize] = -20;
	      if(sigtype!=0)
		{
		  Gmu1pt[gsize] = GenInfo_pt[GenInfo_da1[GenInfo_da1[j]]];
		  Gmu1eta[gsize] = GenInfo_eta[GenInfo_da1[GenInfo_da1[j]]];
		  Gmu1phi[gsize] = GenInfo_phi[GenInfo_da1[GenInfo_da1[j]]];
		  Gmu1p[gsize] = Gmu1pt[gsize]*cosh(Gmu1eta[gsize]);
		  Gmu2pt[gsize] = GenInfo_pt[GenInfo_da2[GenInfo_da1[j]]];
		  Gmu2eta[gsize] = GenInfo_eta[GenInfo_da2[GenInfo_da1[j]]];
		  Gmu2phi[gsize] = GenInfo_phi[GenInfo_da2[GenInfo_da1[j]]];
		  Gmu2p[gsize] = Gmu2pt[gsize]*cosh(Gmu2eta[gsize]);
		  if(sigtype==1||sigtype==2)
		    {
		      Gtk1pt[gsize] = GenInfo_pt[GenInfo_da2[j]];
		      Gtk1eta[gsize] = GenInfo_eta[GenInfo_da2[j]];
		      Gtk1phi[gsize] = GenInfo_phi[GenInfo_da2[j]];
		    }
		  else
		    {
		      Gtk1pt[gsize] = GenInfo_pt[GenInfo_da1[GenInfo_da2[j]]];
		      Gtk1eta[gsize] = GenInfo_eta[GenInfo_da1[GenInfo_da2[j]]];
		      Gtk1phi[gsize] = GenInfo_phi[GenInfo_da1[GenInfo_da2[j]]];
		      Gtk2pt[gsize] = GenInfo_pt[GenInfo_da2[GenInfo_da2[j]]];
		      Gtk2eta[gsize] = GenInfo_eta[GenInfo_da2[GenInfo_da2[j]]];
		      Gtk2phi[gsize] = GenInfo_phi[GenInfo_da2[GenInfo_da2[j]]];
		    }
		}
	      gsize++;
	    }
	  ntGen->Fill();
	}
    }

  outf->Write();
  outf->Close();
  return 1;
}
开发者ID:HyunchulKim,项目名称:BntupleRunII,代码行数:101,代码来源:loop.C

示例3: Loop

void oniaEff::Loop(const char* fname, bool ispbpb, bool isPsip)
{
//   In a ROOT session, you can do:
//      root> .L oniaEff.C
//      root> oniaEff t
//      root> t.GetEntry(12); // Fill t data members with entry number 12
//      root> t.Show();       // Show values of entry 12
//      root> t.Show(16);     // Read and show values of entry 16
//      root> t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
   fChain->SetBranchStatus("*",0);  // disable all branches
   fChain->SetBranchStatus("Centrality",1);  
   fChain->SetBranchStatus("HLTriggers",1);  
   fChain->SetBranchStatus("Reco_QQ_*",1);  
   fChain->SetBranchStatus("Gen_QQ_*",1);  
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch
   if (fChain == 0) return;

   Long64_t nentries = fChain->GetEntriesFast();
   TFile *f = new TFile(fname, "RECREATE");
   f->cd();

   // define the histos
   // we need 2 numerators (with and without ctau3d cut), 1 denominator. Let's make them in arrays
   TH1F *hnum_centmid = new TH1F("hnum_centmid","hnum_centmid",nbins_centmid,bins_centmid);
   TH1F *hnum_centfwd = new TH1F("hnum_centfwd","hnum_centfwd",nbins_centfwd,bins_centfwd);
   TH1F *hnum_ptmid = new TH1F("hnum_ptmid","hnum_ptmid",nbins_ptmid,bins_ptmid);
   TH1F *hnum_ptfwd = new TH1F("hnum_ptfwd","hnum_ptfwd",nbins_ptfwd,bins_ptfwd);
   TH1F *hnumcut_centmid = new TH1F("hnumcut_centmid","hnumcut_centmid",nbins_centmid,bins_centmid);
   TH1F *hnumcut_centfwd = new TH1F("hnumcut_centfwd","hnumcut_centfwd",nbins_centfwd,bins_centfwd);
   TH1F *hnumcut_ptmid = new TH1F("hnumcut_ptmid","hnumcut_ptmid",nbins_ptmid,bins_ptmid);
   TH1F *hnumcut_ptfwd = new TH1F("hnumcut_ptfwd","hnumcut_ptfwd",nbins_ptfwd,bins_ptfwd);
   TH1F *hnumptdepcut_centmid = new TH1F("hnumptdepcut_centmid","hnumptdepcut_centmid",nbins_centmid,bins_centmid);
   TH1F *hnumptdepcut_centfwd = new TH1F("hnumptdepcut_centfwd","hnumptdepcut_centfwd",nbins_centfwd,bins_centfwd);
   TH1F *hnumptdepcut_ptmid = new TH1F("hnumptdepcut_ptmid","hnumptdepcut_ptmid",nbins_ptmid,bins_ptmid);
   TH1F *hnumptdepcut_ptfwd = new TH1F("hnumptdepcut_ptfwd","hnumptdepcut_ptfwd",nbins_ptfwd,bins_ptfwd);
   TH1F *hden_centmid = new TH1F("hden_centmid","hden_centmid",nbins_centmid,bins_centmid);
   TH1F *hden_centfwd = new TH1F("hden_centfwd","hden_centfwd",nbins_centfwd,bins_centfwd);
   TH1F *hden_ptmid = new TH1F("hden_ptmid","hden_ptmid",nbins_ptmid,bins_ptmid);
   TH1F *hden_ptfwd = new TH1F("hden_ptfwd","hden_ptfwd",nbins_ptfwd,bins_ptfwd);

   // also store more finely binned histos
   TH1F *hcentfine = new TH1F("hcentfine","hcentfine",200,0,200);
   TH2F *hnum2d = new TH2F("hnum2d","hnum2d",48,0,2.4,150,0,30);
   TH2F *hden2d = new TH2F("hden2d","hden2d",48,0,2.4,150,0,30);

   //nentries = 200;
   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      // if (ientry>1000) break;

      // skip the event if Gen_QQ_size != 1 for now
      b_Gen_QQ_size->GetEntry(ientry); //read only this branch
      if (Gen_QQ_size !=1) continue;

      // apply gen acceptance cuts
      b_Gen_QQ_mupl_4mom->GetEntry(ientry); //read only this branch
      b_Gen_QQ_mumi_4mom->GetEntry(ientry); //read only this branch
      TLorentzVector *tlvgenpl = (TLorentzVector*) Gen_QQ_mupl_4mom->At(0);
      TLorentzVector *tlvgenmi = (TLorentzVector*) Gen_QQ_mumi_4mom->At(0);
      //if (!isGlobalMuonInAccept2015(tlvgenpl) || !isGlobalMuonInAccept2015(tlvgenmi)) continue;
      bool isgenok = true;
      if (!isGlobalMuonInAccept2015(tlvgenpl) || !isGlobalMuonInAccept2015(tlvgenmi)) isgenok=false;
      
      // fill denominators
      b_Gen_QQ_4mom->GetEntry(ientry); //read only this branch
      TLorentzVector *tlvgenqq = (TLorentzVector*) Gen_QQ_4mom->At(0);
      double genpt = tlvgenqq->Pt();
      b_Centrality->GetEntry(ientry); //read only this branch

      // check that the dimuon is inside the analysis bins
      bool gen_inbin = false;
      if (fabs(tlvgenqq->Rapidity())<1.6 && tlvgenqq->Pt()>6.5 && tlvgenqq->Pt()<30.) gen_inbin = true;
      if (fabs(tlvgenqq->Rapidity())>1.6 && fabs(tlvgenqq->Rapidity())<2.4 && tlvgenqq->Pt()>3. && tlvgenqq->Pt()<30.) gen_inbin = true;
      //if (!gen_inbin) continue;
      if (!gen_inbin) isgenok=false;
      
      double weight = ispbpb ? fChain->GetWeight()*findNcoll(Centrality) : 1.;
      double weight_DataMC = 0.0;
      if(fabs(tlvgenqq->Rapidity())<1.6){
	weight_DataMC = 1.659 - 0.0844*genpt + 0.0019*genpt*genpt;
	if(ispbpb && !isPsip) weight_DataMC = 1.88 - 0.120*genpt + 0.0033*genpt*genpt; 
	if(!ispbpb && isPsip) weight_DataMC = 1.059 + 0.0001*genpt - 0.0005*genpt*genpt;
	if(ispbpb && isPsip) weight_DataMC = 0.254 + 0.109*genpt - 0.0033*genpt*genpt;
	
      }
//.........这里部分代码省略.........
开发者ID:CMS-HIN-dilepton,项目名称:DimuonCADIs,代码行数:101,代码来源:oniaEff.C

示例4: Loop

void test_allEvt::Loop()
{
  std::map<std::string, TH1D*> h_1d;
  std::map<std::string, TH2D*> h_2d;

  GoToHXframe = new HTransformToHelicityFrame();
  GoToCSframe = new HTransformToCS();
  
  TFile *fin_ewk_syst = new TFile("fout_invariant_mass_qed.root");
  TH1D *hvpt_ewk_syst = (TH1D*) fin_ewk_syst->Get("vpt_ewk_p8std_div_ewk_p8bad");
  TH1D *hmass_ewk_syst = (TH1D*) fin_ewk_syst->Get("ewk_p8std_div_ewk_p8bad");
  TH1D *hmass_ewk_syst_photos = (TH1D*) fin_ewk_syst->Get("ewk_p8photos_div_ewk_p8bad");
  
  // TH1D *hmass=new TH1D("hmass","hmass",150,50,200);
  // hmass->Sumw2();
  
  // TH1D *hvpt=new TH1D("hvpt","hvpt",300,0,300);
  // hvpt->Sumw2();
  
  if (fChain == 0) return;

  Long64_t nentries = fChain->GetEntriesFast();

  Long64_t nbytes = 0, nb = 0;
  // for (Long64_t jentry=0; jentry<nentries;jentry++) {
  for (Long64_t jentry=0; jentry<10e6;jentry++) {
    // for (Long64_t jentry=0; jentry<10e3;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    // if (Cut(ientry) < 0) continue;
    
    if(jentry%100000==0) std::cout << "jentry= "<< jentry << "/" << nentries << std::endl;
    
    if(ZGen_PostFSR_mass>50
        && MuPosGenStatus1_pt>25
        && MuNegGenStatus1_pt>25
        && TMath::Abs(MuPosGenStatus1_eta)<2.5
        && TMath::Abs(MuNegGenStatus1_eta)<2.5
        ){
          
      muPosGen_status1.SetPtEtaPhiM(MuPosGenStatus1_pt,MuPosGenStatus1_eta,MuPosGenStatus1_phi,MuPosGenStatus1_mass);
      muNegGen_status1.SetPtEtaPhiM(MuNegGenStatus1_pt,MuNegGenStatus1_eta,MuNegGenStatus1_phi,MuNegGenStatus1_mass);
      ZGen_status1 = muPosGen_status1 + muNegGen_status1;
      // mtv = sqrt(2*ptl1*ptl2*(1d0-cos(delphi)))
      double dphi = MuPosGenStatus1_phi - MuNegGenStatus1_phi; 
      if ( dphi > TMath::Pi() ) {
        dphi -= 2.0*TMath::Pi();
      } else if ( dphi <= -TMath::Pi() ) {
        dphi += 2.0*TMath::Pi();
      }
      double ZGen_status1_mT = TMath::Sqrt(2*MuPosGenStatus1_pt*MuNegGenStatus1_pt*(1-TMath::Cos(dphi)));
      // cout << "ZGen_status1_mT= " << ZGen_status1_mT << " ZGen_status1.Mt()= " << ZGen_status1.Mt() << endl;
      
      double vpt_rew_factor = hvpt_ewk_syst->Interpolate(ZGen_status1.Pt());
      double vmass_std_rew_factor = hmass_ewk_syst->Interpolate(ZGen_status1.M());
      double vmass_photos_rew_factor = hmass_ewk_syst_photos->Interpolate(ZGen_status1.M());
          
      // hmass->Fill(ZGen_PostFSR_mass);
      common_stuff::plot1D("hmass", ZGen_PostFSR_mass, 1, h_1d, 150,50,200 );
      common_stuff::plot1D("hNegleptonpt", muNegGen_status1.Pt(), 1, h_1d, 160,25,65 );
      common_stuff::plot1D("hPosleptonpt", muPosGen_status1.Pt(), 1, h_1d, 160,25,65 );
      common_stuff::plot1D("hvpt", ZGen_status1.Pt(), 1, h_1d, 300,0,300 );
      common_stuff::plot1D("hrapidity", ZGen_status1.Rapidity(), 1, h_1d, 25,-2.5,2.5 );
      common_stuff::plot1D("hNegleptoneta", muNegGen_status1.Eta(), 1, h_1d, 25,-2.5,2.5 );
      common_stuff::plot1D("hPosleptoneta", muPosGen_status1.Eta(), 1, h_1d, 25,-2.5,2.5 );
      // common_stuff::plot1D("hmT", ZGen_status1.Mt(), 1, h_1d, 100,50,100 );
      common_stuff::plot1D("hmT", ZGen_status1_mT, 1, h_1d, 100,50,100 );
      
      common_stuff::plot1D("hmass_VpTRewTo_ewk_p8std", ZGen_PostFSR_mass, vpt_rew_factor, h_1d, 150,50,200 );
      common_stuff::plot1D("hNegleptonpt_VpTRewTo_ewk_p8std", muNegGen_status1.Pt(), vpt_rew_factor, h_1d, 160,25,65 );
      common_stuff::plot1D("hPosleptonpt_VpTRewTo_ewk_p8std", muPosGen_status1.Pt(), vpt_rew_factor, h_1d, 160,25,65 );
      common_stuff::plot1D("hvpt_VpTRewTo_ewk_p8std", ZGen_status1.Pt(), vpt_rew_factor, h_1d, 300,0,300 );
      common_stuff::plot1D("hrapidity_VpTRewTo_ewk_p8std", ZGen_status1.Rapidity(), vpt_rew_factor, h_1d, 25,-2.5,2.5 );
      common_stuff::plot1D("hNegleptoneta_VpTRewTo_ewk_p8std", muNegGen_status1.Eta(), vpt_rew_factor, h_1d, 25,-2.5,2.5 );
      common_stuff::plot1D("hPosleptoneta_VpTRewTo_ewk_p8std", muPosGen_status1.Eta(), vpt_rew_factor, h_1d, 25,-2.5,2.5 );
      // common_stuff::plot1D("hmT_VpTRewTo_ewk_p8std", ZGen_status1.Mt(), vpt_rew_factor, h_1d, 100,50,100 );
      common_stuff::plot1D("hmT_VpTRewTo_ewk_p8std", ZGen_status1_mT, vpt_rew_factor, h_1d, 100,50,100 );
      
      common_stuff::plot1D("hmass_VpTRewTo_ewk_p8std_massStd", ZGen_PostFSR_mass, vpt_rew_factor*vmass_std_rew_factor, h_1d, 150,50,200 );
      common_stuff::plot1D("hNegleptonpt_VpTRewTo_ewk_p8std_massStd", muNegGen_status1.Pt(), vpt_rew_factor*vmass_std_rew_factor, h_1d, 160,25,65 );
      common_stuff::plot1D("hPosleptonpt_VpTRewTo_ewk_p8std_massStd", muPosGen_status1.Pt(), vpt_rew_factor*vmass_std_rew_factor, h_1d, 160,25,65 );
      common_stuff::plot1D("hvpt_VpTRewTo_ewk_p8std_massStd", ZGen_status1.Pt(), vpt_rew_factor*vmass_std_rew_factor, h_1d, 300,0,300 );
      common_stuff::plot1D("hrapidity_VpTRewTo_ewk_p8std_massStd", ZGen_status1.Rapidity(), vpt_rew_factor*vmass_std_rew_factor, h_1d, 25,-2.5,2.5 );
      common_stuff::plot1D("hNegleptoneta_VpTRewTo_ewk_p8std_massStd", muNegGen_status1.Eta(), vpt_rew_factor*vmass_std_rew_factor, h_1d, 25,-2.5,2.5 );
      common_stuff::plot1D("hPosleptoneta_VpTRewTo_ewk_p8std_massStd", muPosGen_status1.Eta(), vpt_rew_factor*vmass_std_rew_factor, h_1d, 25,-2.5,2.5 );
      // common_stuff::plot1D("hmT_VpTRewTo_ewk_p8std_massStd", ZGen_status1.Mt(), vpt_rew_factor*vmass_std_rew_factor, h_1d, 100,50,100 );
      common_stuff::plot1D("hmT_VpTRewTo_ewk_p8std_massStd", ZGen_status1_mT, vpt_rew_factor*vmass_std_rew_factor, h_1d, 100,50,100 );
      
      common_stuff::plot1D("hmass_VpTRewTo_ewk_p8std_massPhotos", ZGen_PostFSR_mass, vpt_rew_factor*vmass_photos_rew_factor, h_1d, 150,50,200 );
      common_stuff::plot1D("hNegleptonpt_VpTRewTo_ewk_p8std_massPhotos", muNegGen_status1.Pt(), vpt_rew_factor*vmass_photos_rew_factor, h_1d, 160,25,65 );
      common_stuff::plot1D("hPosleptonpt_VpTRewTo_ewk_p8std_massPhotos", muPosGen_status1.Pt(), vpt_rew_factor*vmass_photos_rew_factor, h_1d, 160,25,65 );
      common_stuff::plot1D("hvpt_VpTRewTo_ewk_p8std_massPhotos", ZGen_status1.Pt(), vpt_rew_factor*vmass_photos_rew_factor, h_1d, 300,0,300 );
      common_stuff::plot1D("hrapidity_VpTRewTo_ewk_p8std_massPhotos", ZGen_status1.Rapidity(), vpt_rew_factor*vmass_photos_rew_factor, h_1d, 25,-2.5,2.5 );
      common_stuff::plot1D("hNegleptoneta_VpTRewTo_ewk_p8std_massPhotos", muNegGen_status1.Eta(), vpt_rew_factor*vmass_photos_rew_factor, h_1d, 25,-2.5,2.5 );
      common_stuff::plot1D("hPosleptoneta_VpTRewTo_ewk_p8std_massPhotos", muPosGen_status1.Eta(), vpt_rew_factor*vmass_photos_rew_factor, h_1d, 25,-2.5,2.5 );
      // common_stuff::plot1D("hmT_VpTRewTo_ewk_p8std_massPhotos", ZGen_status1.Mt(), vpt_rew_factor*vmass_photos_rew_factor, h_1d, 100,50,100 );
      common_stuff::plot1D("hmT_VpTRewTo_ewk_p8std_massPhotos", ZGen_status1_mT, vpt_rew_factor*vmass_photos_rew_factor, h_1d, 100,50,100 );
    
    }
//.........这里部分代码省略.........
开发者ID:12345ieee,项目名称:cmg-cmssw,代码行数:101,代码来源:test_allEvt.C

示例5: fill

void fill(int const kf, TLorentzVector* b, double weight, TLorentzVector const& p1Mom, TLorentzVector const& p2Mom, TVector3 v00)
{
   int const centrality = floor(nCent * gRandom->Rndm());

   TVector3 const vertex = getVertex(centrality);
   // smear primary vertex
   // float const sigmaVertex = sigmaVertexCent[cent];
   // TVector3 const vertex(gRandom->Gaus(0, sigmaVertex), gRandom->Gaus(0, sigmaVertex), gRandom->Gaus(0, sigmaVertex));

   v00 += vertex;

   // smear momentum
   TLorentzVector const p1RMom = smearMom(0, p1Mom);
   TLorentzVector const p2RMom = smearMom(0, p2Mom);

   // smear position
   TVector3 const p1RPos = smearPosData(0, vertex.z(), centrality, p1RMom, v00);
   TVector3 const p2RPos = smearPosData(0, vertex.z(), centrality, p2RMom, v00);
   // TVector3 const kRPos = smearPos(kMom, kRMom, v00);
   // TVector3 const pRPos = smearPos(pMom, pRMom, v00);

   // reconstruct
   TLorentzVector const rMom = p1RMom + p2RMom;
   float const p1Dca = dca(p1Mom.Vect(), v00, vertex);
   float const p2Dca = dca(p2Mom.Vect(), v00, vertex);
   float const p1RDca = dca(p1RMom.Vect(), p1RPos, vertex);
   float const p2RDca = dca(p2RMom.Vect(), p2RPos, vertex);

   TVector3 v0;
   float const dca12 = dca1To2(p1RMom.Vect(), p1RPos, p2RMom.Vect(), p2RPos, v0);
   float const decayLength = (v0 - vertex).Mag();
   float const dcaD0ToPv = dca(rMom.Vect(), v0, vertex);
   float const cosTheta = (v0 - vertex).Unit().Dot(rMom.Vect().Unit());
   float const angle12 = p1RMom.Vect().Angle(p2RMom.Vect());

   TLorentzVector p1RMomRest = p1RMom;
   TVector3 beta;
   beta.SetMagThetaPhi(rMom.Beta(), rMom.Theta(), rMom.Phi());
   p1RMomRest.Boost(-beta);
   float const cosThetaStar = rMom.Vect().Unit().Dot(p1RMomRest.Vect().Unit());

   // save
   float arr[100];
   int iArr = 0;
   arr[iArr++] = centrality;
   arr[iArr++] = vertex.X();
   arr[iArr++] = vertex.Y();
   arr[iArr++] = vertex.Z();

   arr[iArr++] = kf;
   arr[iArr++] = b->M();
   arr[iArr++] = b->Perp();
   arr[iArr++] = b->PseudoRapidity();
   arr[iArr++] = b->Rapidity();
   arr[iArr++] = b->Phi();
   arr[iArr++] = v00.X();
   arr[iArr++] = v00.Y();
   arr[iArr++] = v00.Z();

   arr[iArr++] = rMom.M();
   arr[iArr++] = rMom.Perp();
   arr[iArr++] = rMom.PseudoRapidity();
   arr[iArr++] = rMom.Rapidity();
   arr[iArr++] = rMom.Phi();
   arr[iArr++] = v0.X();
   arr[iArr++] = v0.Y();
   arr[iArr++] = v0.Z();

   arr[iArr++] = dca12;
   arr[iArr++] = decayLength;
   arr[iArr++] = dcaD0ToPv;
   arr[iArr++] = cosTheta;
   arr[iArr++] = angle12;
   arr[iArr++] = cosThetaStar;

   arr[iArr++] = p1Mom.M();
   arr[iArr++] = p1Mom.Perp();
   arr[iArr++] = p1Mom.PseudoRapidity();
   arr[iArr++] = p1Mom.Rapidity();
   arr[iArr++] = p1Mom.Phi();
   arr[iArr++] = p1Dca;

   arr[iArr++] = p1RMom.M();
   arr[iArr++] = p1RMom.Perp();
   arr[iArr++] = p1RMom.PseudoRapidity();
   arr[iArr++] = p1RMom.Rapidity();
   arr[iArr++] = p1RMom.Phi();
   arr[iArr++] = p1RPos.X();
   arr[iArr++] = p1RPos.Y();
   arr[iArr++] = p1RPos.Z();
   arr[iArr++] = p1RDca;
   arr[iArr++] = tpcReconstructed(0,1,centrality,p1RMom);

   arr[iArr++] = p2Mom.M();
   arr[iArr++] = p2Mom.Perp();
   arr[iArr++] = p2Mom.PseudoRapidity();
   arr[iArr++] = p2Mom.Rapidity();
   arr[iArr++] = p2Mom.Phi();
   arr[iArr++] = p2Dca;

//.........这里部分代码省略.........
开发者ID:amcw7777,项目名称:auau200GeVRun14Ana,代码行数:101,代码来源:toyMcEffKS.C

示例6: fill

void fill(int const kf, TLorentzVector* b, double const weight, TClonesArray& daughters)
{
   TLorentzVector kMom;
   TLorentzVector p1Mom;
   TLorentzVector p2Mom;

   int nTrk = daughters.GetEntriesFast();
   for (int iTrk = 0; iTrk < nTrk; ++iTrk)
   {
      TParticle* ptl0 = (TParticle*)daughters.At(iTrk);

      switch(abs(ptl0->GetPdgCode()))
      {
        case 321:
          ptl0->Momentum(kMom);
          break;
        case 211:
          if(!p1Mom.P()) ptl0->Momentum(p1Mom);
          else ptl0->Momentum(p2Mom);
          break;
        default:
          break;
      }
   }
   daughters.Clear();

   // smear and get total momentum
   TLorentzVector kRMom = smearMom(kMom,fKaonMomResolution);
   TLorentzVector p1RMom = smearMom(p1Mom,fPionMomResolution);
   TLorentzVector p2RMom = smearMom(p2Mom,fPionMomResolution);
   TLorentzVector rMom = kRMom + p1RMom + p2RMom;

   // save
   float arr[100];
   int iArr = 0;
   arr[iArr++] = kf;
   arr[iArr++] = weight;
   arr[iArr++] = b->M();
   arr[iArr++] = b->Perp();
   arr[iArr++] = b->PseudoRapidity();
   arr[iArr++] = b->Rapidity();
   arr[iArr++] = b->Phi();

   arr[iArr++] = rMom.M();
   arr[iArr++] = rMom.Perp();
   arr[iArr++] = rMom.PseudoRapidity();
   arr[iArr++] = rMom.Rapidity();
   arr[iArr++] = rMom.Phi();

   arr[iArr++] = kMom.M();
   arr[iArr++] = kMom.Perp();
   arr[iArr++] = kMom.PseudoRapidity();
   arr[iArr++] = kMom.Rapidity();
   arr[iArr++] = kMom.Phi();

   arr[iArr++] = kRMom.M();
   arr[iArr++] = kRMom.Perp();
   arr[iArr++] = kRMom.PseudoRapidity();
   arr[iArr++] = kRMom.Rapidity();
   arr[iArr++] = kRMom.Phi();

   arr[iArr++] = p1Mom.M();
   arr[iArr++] = p1Mom.Perp();
   arr[iArr++] = p1Mom.PseudoRapidity();
   arr[iArr++] = p1Mom.Rapidity();
   arr[iArr++] = p1Mom.Phi();

   arr[iArr++] = p1RMom.M();
   arr[iArr++] = p1RMom.Perp();
   arr[iArr++] = p1RMom.PseudoRapidity();
   arr[iArr++] = p1RMom.Rapidity();
   arr[iArr++] = p1RMom.Phi();

   arr[iArr++] = p2Mom.M();
   arr[iArr++] = p2Mom.Perp();
   arr[iArr++] = p2Mom.PseudoRapidity();
   arr[iArr++] = p2Mom.Rapidity();
   arr[iArr++] = p2Mom.Phi();

   arr[iArr++] = p2RMom.M();
   arr[iArr++] = p2RMom.Perp();
   arr[iArr++] = p2RMom.PseudoRapidity();
   arr[iArr++] = p2RMom.Rapidity();
   arr[iArr++] = p2RMom.Phi();
   
   nt->Fill(arr);
}
开发者ID:GuannanXie,项目名称:auau200GeVRun14,代码行数:87,代码来源:pythia_decay.C

示例7: RunPidGetterQAEff

void RunPidGetterQAEff()
{
        
    TString PidFrameworkDir = "/lustre/nyx/cbm/users/klochkov/soft/PidFramework/";
    gSystem->Load( PidFrameworkDir + "build/libPid");  
    
    gStyle->SetOptStat(0000);

    TFile *f2 = new TFile("pid_0.root");    
    TTree *PidTree = (TTree*) f2->Get("PidTree");

    TofPidGetter *getter = new TofPidGetter();
    TBranch *PidGet = PidTree->GetBranch("TofPidGetter");
    PidGet->SetAddress(&getter);

    PidGet->GetEntry(0);
    Float_t ret[3];
        
    TProfile *hEffP  = new  TProfile ("hEffP", "", 100, 0, 10);
    TProfile *hEffPi = new  TProfile ("hEffPi", "", 100, 0, 6);
    TProfile *hEffK  = new  TProfile ("hEffK", "", 100, 0, 5);

    TProfile *hEffPSigma  = new  TProfile ("hEffPSigma", "", 100, 0, 10);
    TProfile *hEffPiSigma = new  TProfile ("hEffPiSigma", "", 100, 0, 6);
    TProfile *hEffKSigma  = new  TProfile ("hEffKSigma", "", 100, 0, 5);

    TProfile2D *hEffPtYP  = new  TProfile2D ("hEffPtYP", "", 100,   -2.5, 2.5, 100, 0, 4);
    TProfile2D *hEffPtYK  = new  TProfile2D ("hEffPtYK", "", 100,   -2.5, 2.5, 100, 0, 4);
    TProfile2D *hEffPtYPi  = new  TProfile2D ("hEffPtYPi", "", 100, -2.5, 2.5, 100, 0, 4);

    TString InTreeFileName = "/lustre/nyx/cbm/users/dblau/cbm/mc/UrQMD/AuAu/10AGeV/sis100_electron/SC_ON/2016_09_01/tree/11111.root";
    TFile *InFile = new TFile(InTreeFileName);    
    TTree *InTree = (TTree*) InFile->Get("fDataTree");
    DataTreeEvent* DTEvent;
    InTree -> SetBranchAddress("DTEvent",&DTEvent);

    int nevents = 100000;//InTree->GetEntries();
    int outputstep = 100;
    std::cout << "Entries = " << nevents << std::endl;

    for (int j=0;j<nevents;j++)
    {
        if ( (j+1) % outputstep == 0) std::cout << j+1 << "/" << nevents <<  "\r" << std::flush;
        InTree->GetEntry(j);
        
        Int_t Nmc[3] = {100,100,100};
        Int_t Ntof[3] = {0,0,0};
        Int_t PdgCode[3] = {2212, 212, 211};
        Double_t sigmas [3] = {0,0,0};
                
        for (int i=0;i<DTEvent->GetNTracks(); i++)
        {
            
            TLorentzVector v;
            
            DataTreeTrack* track = DTEvent -> GetTrack(i);
            DataTreeMCTrack* mctrack = DTEvent -> GetMCTrack(i);
            Double_t p = mctrack->GetPt() * TMath::CosH( mctrack->GetEta() );
            if (track->GetTOFHitId() < 0)
            {
                if (mctrack->GetPdgId() == 2212  )  
                {
                    v.SetPtEtaPhiM (track->GetPt(0), track->GetEta(0), track->GetPhi(0), 0.9386);
                    hEffP->Fill ( p, 0 );
                    hEffPSigma->Fill ( p, 0 );
                    hEffPtYP -> Fill( v.Rapidity() - 1.52, v.Pt(), 0 );
                }
                
                if (mctrack->GetPdgId() == 321  )  
                {
                    v.SetPtEtaPhiM (track->GetPt(0), track->GetEta(0), track->GetPhi(0), 0.5);
                    hEffK->Fill ( p, 0 );
                    hEffKSigma->Fill ( p,  0 ); 
                    hEffPtYK -> Fill( v.Rapidity() - 1.52, v.Pt(), 0 );
                }
                
                if (mctrack->GetPdgId() == 211  )  
                {
                    v.SetPtEtaPhiM (track->GetPt(0), track->GetEta(0), track->GetPhi(0), 0.14);
                    hEffPi->Fill ( p, 0 );
                    hEffPiSigma->Fill ( p,  0);  
                    hEffPtYPi -> Fill( v.Rapidity() - 1.52, v.Pt(),  0 );
                }
                continue;
            }
//             
            DataTreeTOFHit* toftrack = DTEvent -> GetTOFHit(track->GetTOFHitId());
            p = toftrack->GetP();
            Double_t m2 = toftrack->GetMass2 ();
            
            Bool_t cut = toftrack->GetBeta() > 0.1 && ( track->GetChiSq(0)/track->GetNDF(0) < 3 ) && p > 1.0 ;
            if ( !cut ) continue;
            
            getter->GetBayesProbability (m2, p, ret);
            sigmas[0] = getter->GetSigmaProton (m2, p);
            sigmas[1] = getter->GetSigmaKaon (m2, p);
            sigmas[2] = getter->GetSigmaPion (m2, p);
            
    //             std::cout << "pdg = " << mctrack->GetPdgId() << " p = " << p << " Sp = " << sigmas[0] << " Sk = " << sigmas[1]<< " Spi = " << sigmas[2] << std::endl;
            
//.........这里部分代码省略.........
开发者ID:viktorklochkov,项目名称:Pid,代码行数:101,代码来源:RunPidGetterQAEff.C

示例8: createWorkspace


//.........这里部分代码省略.........

    	const double bordersPT[nPT+1] = {0., 12., 16., 20., 30., 100.};
    	const double bordersRAP[nRAP+1] = {0., 0.6, 2.};
    	const double bordersL[nL+1] = {onia::ctVarMin, -0.05, -0.03, -0.02, -0.015, -0.01, -0.005,  0., 0.005, 0.01, 0.015, 0.02, 0.03, 0.05, 0.1, onia::ctVarMax};

    	TH1D* h_ctauerr_2011[nRAP+1][nPT+1][nL+1];
    	TH1D* h_ctauerr_2012[nRAP+1][nPT+1][nL+1];

    	for(int iRAP = 0; iRAP < nRAP+1; iRAP++){
    		for(int iPT = 0; iPT < nPT+1; iPT++){
    			for(int iL = 0; iL < nL+1; iL++){

    				h_ctauerr_2011[iRAP][iPT][iL] = (TH1D*)infile->Get(Form("h_ctauerr_2011_rap%d_pt%d_l%d",iRAP, iPT, iL));
    				h_ctauerr_2012[iRAP][iPT][iL] = (TH1D*)infile->Get(Form("h_ctauerr_2012_rap%d_pt%d_l%d",iRAP, iPT, iL));

    			}
    		}
    	}
    	cout<<"opened hists"<<endl;

    	/// Finished reading in 2011 data ctauErr-histos
    */


    // loop through events in tree and save them to dataset
    for (int ientries = 0; ientries < entries; ientries++) {
        numEntriesTotal++;
        if (ientries%10000==0) std::cout << "event " << ientries << " of " << entries <<  std::endl;

        tree->GetEntry(ientries);

        double M_jpsi =jpsi->M();
        double pt_jpsi =jpsi->Pt();
        double y_jpsi =jpsi->Rapidity();
        double M =chic_rf->M();
        //double M =chic->M()-jpsi->M()+onia::MpsiPDG;
        double y=chic->Rapidity();
        double pt=chic->Pt();

        //if (ientries%3==0){
        //	M_jpsi = gRandom->Uniform(JpsiMass->getMin(), JpsiMass->getMax());
        //}


        //double JpsictErrRand = h_JpsictErr->GetRandom();
        //double JpsictErrRand2 = h_JpsictErr->GetRandom();
        //double JpsictMeanRand=0.;
        ////double pTcorrection=(pt-20.)*0.002;
        //
        ////JpsictErrRand-=pTcorrection;
        //if(JpsictErrRand<0) JpsictErrRand=0.001;
        ////JpsictErrRand2-=pTcorrection;
        //if(JpsictErrRand2<0) JpsictErrRand2=0.001;
        //
        //if (ientries%1000000==0){
        //	double exponent=0.4;
        //	JpsictMeanRand=gRandom->Exp(exponent);
        //}
        //
        //
        //lifetime = gRandom->Gaus(JpsictMeanRand,0.8*JpsictErrRand);
        //lifetimeErr = JpsictErrRand2;
        //if (ientries%3==0){
        //	lifetime = gRandom->Gaus(JpsictMeanRand,1.5*JpsictErrRand);
        //}
        //
开发者ID:knuenz,项目名称:ChicPol,代码行数:67,代码来源:createWorkspace.C

示例9: loopNonpromptBzero


//.........这里部分代码省略.........
  int type,flag;
  int flagEvt=0;  
  int offsetHltTree=0;

  int testevent=0,testcand=0;
  
  for (Long64_t i=0; i<nentries;i++) {
    nbytes += root->GetEntry(i);
    flagEvt=0;
    while (flagEvt==0)
    {
       hlt->GetEntry(i+offsetHltTree);
       //cout <<offsetHltTree<<" "<<Bfr_HLT_Event<<" "<<EvtInfo_EvtNo<<endl;
       if (Bfr_HLT_Event==EvtInfo_EvtNo && Bfr_HLT_Run==EvtInfo_RunNo) flagEvt=1; else offsetHltTree++;
    } 

    if (i%10000==0) cout <<i<<" / "<<nentries<<"   offset HLT:"<<offsetHltTree<<endl;

    int type1size=0,type2size=0,type3size=0,type4size=0,type5size=0,type6size=0,type7size=0;
    float best,best2,temy;
    int bestindex,best2index;

    size=0;
    best=-1;
    bestindex=-1;
    best2=10000.;
    best2index=-1;
    for (int j=0;j<BInfo_size;j++) 
      {
	if(BInfo_type[j]>7) continue;
	if (ifchannel[BInfo_type[j]-1]!=1) continue;
	//skim{{{
	b4Pout->SetXYZM(BInfo_px[j],BInfo_py[j],BInfo_pz[j],BInfo_mass[j]);
	temy = b4Pout->Rapidity();
	if(REAL)
	  {
	    if(!(((EvtInfo_RunNo>=210498&&EvtInfo_RunNo<=211256&&abs(temy+0.465)<1.93)||(EvtInfo_RunNo>=211313&&EvtInfo_RunNo<=211631&&abs(temy-0.465)<1.93)))) continue;
	  }
	else
	  {
	    if(abs(temy+0.465)>=1.93) continue;
	  }
	if(BInfo_mass[j]<5 || BInfo_mass[j]>6) continue;
	if(BInfo_pt[j]<10.) continue;
	//}}}
	if(BInfo_type[j]==1)
	  {
	    //if(TrackInfo_pt[BInfo_rftk1_index[j]]<0.9) continue;
	    fillTree(bP,bVtx,b4P,j,type1size,KAON_MASS,0,REAL);
	    if(chi2cl[type1size]>best)
	      {
		best = chi2cl[type1size];
		bestindex = type1size;
	      }
	    type1size++;
	  }
      }
    if(size>0)
    {
      bestchi2 = bestindex;
      isbestchi2[bestindex] = 1;
    }
    nt0->Fill();

    size=0;
    best=-1;
开发者ID:KiSooLee,项目名称:Bntuple,代码行数:67,代码来源:loopNonpromptBzero.C

示例10: analyzer


//.........这里部分代码省略.........
        continue;

    /////////////////////////////////////////////////
    // Muon Resolution Smearing to match MuScleFit

    if(isSignal) // smear only signal because it has muons from higgs 
    {
      if(reco1GenPostFSR.pt<0.)
        cout << "Muon 1 Post FSR not valid!\n";
      if(reco2GenPostFSR.pt<0.)
        cout << "Muon 2 Post FSR not valid!\n";
      float ptReco1 = -1.;
      float ptReco2 = -1.;
      if(runPeriod == "7TeV")
      {
        ptReco1 = smearPT2011 -> PTsmear(reco1GenPostFSR.pt, reco1GenPostFSR.eta, reco1GenPostFSR.charge, reco1.pt, ISMEAR);
        ptReco2 = smearPT2011 -> PTsmear(reco2GenPostFSR.pt, reco2GenPostFSR.eta, reco2GenPostFSR.charge, reco2.pt, ISMEAR);
      }
      else
      {
        ptReco1 = smearPT -> PTsmear(reco1GenPostFSR.pt, reco1GenPostFSR.eta, reco1GenPostFSR.charge, reco1.pt, ISMEAR);
        ptReco2 = smearPT -> PTsmear(reco2GenPostFSR.pt, reco2GenPostFSR.eta, reco2GenPostFSR.charge, reco2.pt, ISMEAR);
      }
      TLorentzVector reco1Vec;
      TLorentzVector reco2Vec;
      reco1Vec.SetPtEtaPhiM(ptReco1,reco1.eta,reco1.phi,MASS_MUON);
      reco2Vec.SetPtEtaPhiM(ptReco2,reco2.eta,reco2.phi,MASS_MUON);
      TLorentzVector diMuonVec = reco1Vec + reco2Vec;

      reco1.pt = ptReco1;
      reco2.pt = ptReco2;
      recoCandMass = diMuonVec.M();
      recoCandPt = diMuonVec.Pt();
      recoCandY = diMuonVec.Rapidity();
      recoCandPhi = diMuonVec.Phi();
    
      reco1Vec.SetPtEtaPhiM(ptReco1,reco1.eta,reco1.phi,MASS_MUON);
      reco2Vec.SetPtEtaPhiM(ptReco2,reco2.eta,reco2.phi,MASS_MUON);
      diMuonVec = reco1Vec + reco2Vec;
    }

    //////////////////////////////////////////
    // Muon-related cuts

    if (recoCandMass > maxMmm || recoCandMass < minMmm)
        continue;

    bool muon1PassId = (*muonIdFuncPtr)(reco1);
    bool muon2PassId = (*muonIdFuncPtr)(reco2);
    if ( !(muon1PassId && muon2PassId))
        continue;

    bool muon1PassIso = (getPFRelIso(reco1) <= 0.12);
    bool muon2PassIso = (getPFRelIso(reco2) <= 0.12);
    if ( !(muon1PassIso && muon2PassIso))
        continue;

    // Order muons by pt
    if (reco1.pt < reco2.pt)
    {
      _MuonInfo tmpMuon = reco1;
      reco1 = reco2;
      reco1 = tmpMuon;
    }

    // PU reweight
开发者ID:jhugon,项目名称:hmumuToyAnalysis,代码行数:67,代码来源:analyzer.C

示例11: createWorkspace

void createWorkspace(const std::string &infilename, int nState, bool correctCtau, bool drawRapPt2D, bool drawPtCPM2D){
	gROOT->SetStyle("Plain");
	gStyle->SetTitleBorderSize(0);

	// Set some strings
	const std::string workspacename = "ws_masslifetime",
				treename = "selectedData";

	// Get the tree from the data file
	TFile *f = TFile::Open(infilename.c_str());
	TTree *tree = (TTree*)f->Get(treename.c_str());

	// Set branch addresses in tree to be able to import tree to roofit
	TLorentzVector* jpsi = new TLorentzVector;
	tree->SetBranchAddress("JpsiP",&jpsi);
	double CPMval = 0;
	tree->SetBranchAddress("CPM",&CPMval);
	double massErr = 0;
	tree->SetBranchAddress("JpsiMassErr",&massErr);
	double Vprob = 0;
	tree->SetBranchAddress("JpsiVprob",&Vprob);
	double lifetime = 0;
	tree->SetBranchAddress("Jpsict",&lifetime);
	double lifetimeErr = 0;
	tree->SetBranchAddress("JpsictErr",&lifetimeErr);

	// define variables necessary for J/Psi(Psi(2S)) mass,lifetime fit
	RooRealVar* JpsiMass =
		new RooRealVar("JpsiMass", "M [GeV]", onia::massMin, onia::massMax);
	RooRealVar* JpsiMassErr =
		new RooRealVar("JpsiMassErr", "#delta M [GeV]", 0, 5);
	RooRealVar* JpsiRap =
		new RooRealVar("JpsiRap", "y", -onia::rap, onia::rap);
	RooRealVar* JpsiPt =
		new RooRealVar("JpsiPt", "p_{T} [GeV]", 0. ,100.);
	RooRealVar* JpsiCPM =
		new RooRealVar("JpsiCPM", "N_{ch}", 0. ,100.);		
	RooRealVar* Jpsict =
		new RooRealVar("Jpsict", "lifetime [mm]", -1., 2.5);
	RooRealVar* JpsictErr =
		new RooRealVar("JpsictErr", "Error on lifetime [mm]", 0.0001, 1);
	RooRealVar* JpsiVprob =
		new RooRealVar("JpsiVprob", "", 0.01, 1.);

	// Set bins
	Jpsict->setBins(10000,"cache");
	Jpsict->setBins(100);
	JpsiMass->setBins(100);
	JpsictErr->setBins(100);

	// The list of data variables    
	RooArgList dataVars(*JpsiMass,*JpsiMassErr,*JpsiRap,*JpsiPt,*JpsiCPM,*Jpsict,*JpsictErr,*JpsiVprob);

	// construct dataset to contain events
	RooDataSet* fullData = new RooDataSet("fullData","The Full Data From the Input ROOT Trees",dataVars);

	int entries = tree->GetEntries();
	cout << "entries " << entries << endl;

	// loop through events in tree and save them to dataset
	for (int ientries = 0; ientries < entries; ientries++) {
	
		if (ientries%100000==0) std::cout << "event " << ientries << " of " << entries <<  std::endl;

		tree->GetEntry(ientries);

		double M =jpsi->M();
		double y=jpsi->Rapidity();
		double pt=jpsi->Pt();
		double cpm=CPMval;


		if (M > JpsiMass->getMin() && M < JpsiMass->getMax()
				&& massErr > JpsiMassErr->getMin() && massErr < JpsiMassErr->getMax()
				&& pt > JpsiPt->getMin() && pt < JpsiPt->getMax()
				&& cpm > JpsiCPM->getMin() && cpm < JpsiCPM->getMax()
				&& y > JpsiRap->getMin() && y < JpsiRap->getMax()
				&& lifetime > Jpsict->getMin() && lifetime < Jpsict->getMax()
				&& lifetimeErr > JpsictErr->getMin() && lifetimeErr < JpsictErr->getMax()
				&& Vprob > JpsiVprob->getMin() && Vprob < JpsiVprob->getMax()
			 ){

			JpsiPt      ->setVal(pt); 
			JpsiCPM		->setVal(cpm);
			JpsiRap     ->setVal(y); 
			JpsiMass    ->setVal(M);
			JpsiMassErr ->setVal(massErr);
			JpsiVprob   ->setVal(Vprob);

			//cout<<"before lifetime correction \n"
			//	<<"Jpsict: "<<lifetime<<" JpsictErr: "<<lifetimeErr<<endl;

			if(correctCtau){
				lifetime    = lifetime    * onia::MpsiPDG / M ;
				lifetimeErr = lifetimeErr * onia::MpsiPDG / M ;
				Jpsict    ->setVal(lifetime);
				JpsictErr ->setVal(lifetimeErr);
				//cout<<"MpsiPDG: "<<onia::MpsiPDG<<endl;
				//cout<<"after lifetime correction \n"
				//	<<"Jpsict: "<<lifetime<<" JpsictErr: "<<lifetimeErr<<endl;
//.........这里部分代码省略.........
开发者ID:cferraio,项目名称:JPsi_Nch_Polarization,代码行数:101,代码来源:createWorkspace.C

示例12: Loop

//==============================================
void PolMC::Loop(Int_t selDimuType)
{
  if (fChain == 0) return;

  Long64_t nentries = fChain->GetEntries();
  Long64_t countGenEvent = 0;
  Long64_t nb = 0;
  printf("number of entries = %d\n", (Int_t) nentries);

  //loop over the events
  for (Long64_t jentry=0; jentry<nentries;jentry++) {

    if(jentry % 100000 == 0) printf("event %d\n", (Int_t) jentry);

    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);

    //protection against dummy events:
    if(muPosPx_Gen < -9900.)
      continue;

    hGen_StatEv->Fill(0.5);//count all events
    //do NOT select on the dimuon type (only a RECO variable)
    hGen_StatEv->Fill(1.5);//count all events

    Double_t enMuPos = sqrt(muPosPx_Gen*muPosPx_Gen + muPosPy_Gen*muPosPy_Gen + muPosPz_Gen*muPosPz_Gen + jpsi::muMass*jpsi::muMass);
    Double_t enMuNeg = sqrt(muNegPx_Gen*muNegPx_Gen + muNegPy_Gen*muNegPy_Gen + muNegPz_Gen*muNegPz_Gen + jpsi::muMass*jpsi::muMass);
    TLorentzVector *muPos = new TLorentzVector();
    TLorentzVector *muNeg = new TLorentzVector();
    muPos->SetPxPyPzE(muPosPx_Gen, muPosPy_Gen, muPosPz_Gen, enMuPos);
    muNeg->SetPxPyPzE(muNegPx_Gen, muNegPy_Gen, muNegPz_Gen, enMuNeg);
    Double_t etaMuPos = muPos->PseudoRapidity();
    Double_t etaMuNeg = muNeg->PseudoRapidity();
    Double_t pTMuPos = muPos->Pt();
    Double_t pTMuNeg = muNeg->Pt();

    //take muons only within a certain eta range
    if(TMath::Abs(etaMuPos) > jpsi::etaPS || TMath::Abs(etaMuNeg) > jpsi::etaPS){
      // printf("eta(pos. muon) = %f, eta(neg. muon) = %f\n", etaMuPos, etaMuNeg);
      continue;
    }
    hGen_StatEv->Fill(2.5);//count all events

    if(pTMuPos < jpsi::pTMuMin && pTMuNeg < jpsi::pTMuMin){
      // printf("pT(pos. muon) = %f, pT(neg. muon) = %f\n", pTMuPos, pTMuNeg);
      continue;
    }
    hGen_StatEv->Fill(3.5);

    //test according to Gavin's proposal:
    //if any of the two muons is within 1.4 < eta < 1.6 AND
    //the two muons are close in eta (deltaEta < 0.2)
    //reject the dimuon (no matter whether it is "Seagull" or 
    //"Cowboy"):
//     if(TMath::Abs(etaMuPos - etaMuNeg) < 0.2 &&
//        ((etaMuPos > 1.4 && etaMuPos < 1.6) || (etaMuNeg > 1.4 && etaMuNeg < 1.6))){
//       printf("rejecting the event!\n");
//       continue;
//     }

    //build the invariant mass, pt, ... of the two muons
    TLorentzVector *onia = new TLorentzVector();
    *onia = *(muPos) + *(muNeg);

    Double_t onia_mass = onia->M();
    Double_t onia_pt = onia->Pt();
    Double_t onia_P = onia->P();
    Double_t onia_eta = onia->PseudoRapidity();
    Double_t onia_rap = onia->Rapidity();
    Double_t onia_phi = onia->Phi();
    Double_t onia_mT = sqrt(onia_mass*onia_mass + onia_pt*onia_pt);

    // Int_t rapIndex = -1;
    // for(int iRap = 0; iRap < 2*jpsi::kNbRapBins; iRap++){
    //   if(onia_rap > jpsi::rapRange[iRap] && onia_rap < jpsi::rapRange[iRap+1]){
    // 	rapIndex = iRap+1;
    // 	break;
    //   }
    // }
    Int_t rapForPTIndex = -1;
    for(int iRap = 0; iRap < jpsi::kNbRapForPTBins; iRap++){
      if(TMath::Abs(onia_rap) > jpsi::rapForPTRange[iRap] && 
	 TMath::Abs(onia_rap) < jpsi::rapForPTRange[iRap+1]){
	 rapForPTIndex = iRap+1;
	break;
      }
    }
    Int_t pTIndex = -1;
    for(int iPT = 0; iPT < jpsi::kNbPTBins[rapForPTIndex]; iPT++){
      if(onia_pt > jpsi::pTRange[rapForPTIndex][iPT] && onia_pt < jpsi::pTRange[rapForPTIndex][iPT+1]){
	pTIndex = iPT+1;
	break;
      }
    }
    Int_t rapIntegratedPTIndex = -1;
    for(int iPT = 0; iPT < jpsi::kNbPTBins[0]; iPT++){
      if(onia_pt > jpsi::pTRange[0][iPT] && onia_pt < jpsi::pTRange[0][iPT+1]){
	rapIntegratedPTIndex = iPT+1;
//.........这里部分代码省略.........
开发者ID:hwoehri,项目名称:UserCode,代码行数:101,代码来源:PolMC.C

示例13: dieleAna


//.........这里部分代码省略.........
        eventmixer_eff[multbin].addVector(vep_eff,2);
        eventmixer_eff[multbin].addVector(vem_eff,3);
        vector<pair<HGeantKine *, HGeantKine* > >& pairsVec_eff_multbin = eventmixer_eff[multbin].getMixedVector();

        for (int imix = 0; imix < 3; ++imix) {
            vector<pair<HGeantKine *, HGeantKine* > > pairsVec;
            TString suffix;
            switch (imix) {
                case 0: 
                    pairsVec = pairsVec_acc;
                    suffix = TString(""); 
                    break;
                case 1: 
                    pairsVec = pairsVec_eff;
                    suffix = TString("_eff_multbin0");
                    break;
                case 2: 
                    pairsVec = pairsVec_eff_multbin;
                    suffix = TString("_eff_multbin")+TString::Itoa(multbin,10);
                    break;
            }
            size = pairsVec.size();

            for (Int_t j = 0; j < size; j ++) {
                pair<HGeantKine*,HGeantKine*>& pair = pairsVec[j];

                kine1 = pair.first;
                kine2 = pair.second;
                TLorentzVector vec1, vec2;
                HParticleTool::getTLorentzVector(kine1,vec1,kine1->getID());
                HParticleTool::getTLorentzVector(kine2,vec2,kine2->getID());

                Float_t mom1 = vec1.Vect().Mag();
                Float_t the1 = vec1.Theta()*TMath::RadToDeg();
                Float_t mom2 = vec2.Vect().Mag();
                Float_t the2 = vec2.Theta()*TMath::RadToDeg();

                TLorentzVector dilep = vec1 + vec2;
                Float_t oAngle = vec1.Angle(vec2.Vect())*TMath::RadToDeg();
                Float_t mass   = dilep.M()/1000;
                Float_t pt     = dilep.Perp();
                Float_t y      = dilep.Rapidity();
                Float_t mbinw  = hM.get("hmassNP")->GetBinWidth(hM.get("hmassNP")->FindBin(mass));

                if (oAngle < 9) continue;

                TString chg = "NP";
                if (kine1->getID() == kine2->getID()) {
                    if (kine1->getID() == 2) chg = "PP";
                    if (kine1->getID() == 3) chg = "NN";
                }

                hM.get( TString("hmass")      +chg+suffix)->Fill(mass,       1./mbinw);
                hM.get( TString("hoAngle")    +chg+suffix)->Fill(oAngle      );
                hM.get( TString("hy")         +chg+suffix)->Fill(y           );
                hM.get( TString("hpt")        +chg+suffix)->Fill(pt          );
                hM.get2(TString("hoAnglemass")+chg+suffix)->Fill(oAngle,mass,1./mbinw);
                hM.get2(TString("hoAnglept")  +chg+suffix)->Fill(oAngle,pt   );
                hM.get2(TString("hmasspt")    +chg+suffix)->Fill(mass,pt,    1./mbinw);
                hM.get2(TString("hoAngley")   +chg+suffix)->Fill(oAngle,y    );
                hM.get2(TString("hmassy")     +chg+suffix)->Fill(mass,y,     1./mbinw);
                hM.get2(TString("hpty")       +chg+suffix)->Fill(pt,y        );
                hM.get2(TString("hth1th2")    +chg+suffix)->Fill(the1,the2   );
                hM.get2(TString("hp1p2")      +chg+suffix)->Fill(mom1,mom2   );
            }
        }
//#define DELETE_MIX
#ifdef DELETE_MIX
        vector <HGeantKine *>* toDel = eventmixer.getObjectsToDelete();
        for (unsigned int ii = 0; ii < toDel->size(); ++ii) {
            delete toDel->at(ii);
        }
        toDel->clear();
        delete toDel;
        vector <HGeantKine *>* toDel_eff = eventmixer_eff[0].getObjectsToDelete();
        for (unsigned int ii = 0; ii < toDel_eff->size(); ++ii) {
            delete toDel_eff->at(ii);
        }
        toDel_eff->clear();
        delete toDel_eff;
        vector <HGeantKine *>* toDel_eff_multbin = eventmixer_eff[multbin].getObjectsToDelete();
        for (unsigned int ii = 0; ii < toDel_eff_multbin->size(); ++ii) {
            delete toDel_eff_multbin->at(ii);
        }
        toDel_eff_multbin->clear();
        delete toDel_eff_multbin;
#endif

    } // end event loop

    timer.Stop();

    hM.getFile()->cd();
    TMacro m1(__DIELEANA_FILE__);
    m1.Write();
    hM.writeHists("nomap");

    cout<<"####################################################"<<endl;
    return 0;
}
开发者ID:szymonharabasz,项目名称:hades-programs,代码行数:101,代码来源:dieleAna_em_kine.C

示例14: MakeTree_2015


//.........这里部分代码省略.........

  TClonesArray *Reco_QQ_vtx = 0;;
  TClonesArray *Reco_QQ_4mom = 0;
  TClonesArray *Reco_QQ_mupl_4mom = 0;
  TClonesArray *Reco_QQ_mumi_4mom = 0;
  TClonesArray *Reco_mu_4mom = 0;

  Int_t Reco_mu_size;
  Int_t Reco_mu_type[NMAX];
  Int_t Reco_mu_charge[NMAX];

  int nPlusMu;
  int nMinusMu;
  float muPt[NMAX];
  float muEta[NMAX];
  float muPhi[NMAX];
  //Int_t nReco_QQ;
  Float_t invariantMass;
  Int_t QQtrig;
  Int_t QQsign;
  Int_t QQTrkDr03;
  Int_t QQTrkDr04;
  Int_t QQTrkDr05;

  Int_t QQTrkPt04;
  Int_t QQTrkPt03;
  Int_t QQTrkPt02;

  float weight;
  float weight2;
  Float_t upsPt;
  Float_t upsEta;
  Float_t upsPhi;
  Float_t upsRapidity;	
  Float_t vProb;
  Float_t muPlusPt;
  Float_t muMinusPt;
  Float_t muPlusEta;
  Float_t muMinusEta;
  Float_t muPlusPhi;
  Float_t muMinusPhi;
  Float_t RmuPlusPhi = 0;
  Float_t RmuMinusPhi = 0;

  // centrliaty extra stuff
  Int_t Npix, NpixelTracks,Ntracks;
  Float_t SumET_HF, SumET_HFplus, SumET_HFminus, SumET_HFplusEta4, SumET_HFminusEta4, SumET_EB, SumET_ET, SumET_EE, SumET_EEplus, SumET_EEminus, SumET_ZDC, SumET_ZDCplus, SumET_ZDCminus;

  // track extra stuff
  const int NTRKMAX=100000000;
  int Reco_trk_size=0;
  TClonesArray *Reco_trk_vtx = 0;
  TClonesArray *Reco_trk_4mom = 0;
  float trkPt[NTRKMAX];
  float trkEta[NTRKMAX];
  float trkPhi[NTRKMAX];
  //---------------------------
  TBranch *b_Centrality; //!
  TBranch *b_eventNb;    //!
  TBranch *b_runNb;      //!
  TBranch *b_HLTriggers; //!
  TBranch *b_Reco_QQ_size; //!
  TBranch *b_Reco_QQ_trig; //!
  TBranch *b_Reco_QQ_type; //!
  TBranch *b_Reco_QQ_sign; //!
  TBranch *b_Reco_QQ_VtxProb;  //!
开发者ID:goni,项目名称:HIN-15-001,代码行数:67,代码来源:MakeTree_2015.C

示例15: Loop

void ljpsiresol::Loop(const char* fname, bool isdata, bool ispbpb)
{
//   In a ROOT session, you can do:
//      root> .L ljpsiresol.C
//      root> ljpsiresol t
//      root> t.GetEntry(12); // Fill t data members with entry number 12
//      root> t.Show();       // Show values of entry 12
//      root> t.Show(16);     // Read and show values of entry 16
//      root> t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
   fChain->SetBranchStatus("*",0);  // disable all branches
   fChain->SetBranchStatus("HLTriggers",1);  
   fChain->SetBranchStatus("Centrality",1);  
   fChain->SetBranchStatus("Reco_QQ_*",1);  
   fChain->SetBranchStatus("Gen_QQ_*",!isdata);  
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch
   if (fChain == 0) return;

   Long64_t nentries = fChain->GetEntriesFast();

   TFile *f = new TFile(fname, "RECREATE");
   f->cd();

   // define the histos
   // we need 2 numerators (with and without ctau3d cut), 1 denominator. Let's make them in arrays
   // TH1F *hnum_centmid = new TH1F("hnum_centmid","hnum_centmid",nbins_centmid,bins_centmid);
   vector<TH1F*> hists_mid, hists_fwd, hists_midcent, hists_fwdcent;
   for (int i=0; i<nbins_ptmid; i++) {
      double ptmin = bins_ptmid[i];
      double ptmax = bins_ptmid[i+1];
      TH1F *hist = new TH1F(Form("hmid_pt%.1f-%.1f",ptmin,ptmax),"",70,-0.15,0.2);
      hists_mid.push_back(hist);
   }
   for (int i=0; i<nbins_ptfwd; i++) {
      double ptmin = bins_ptfwd[i];
      double ptmax = bins_ptfwd[i+1];
      TH1F *hist = new TH1F(Form("hfwd_pt%.1f-%.1f",ptmin,ptmax),"",70,-0.15,0.2);
      hists_fwd.push_back(hist);
   }
   for (int i=0; i<nbins_centmid; i++) {
      double centmin = bins_centmid[i];
      double centmax = bins_centmid[i+1];
      TH1F *hist = new TH1F(Form("hmid_cent%.1f-%.1f",centmin,centmax),"",70,-0.15,0.2);
      hists_midcent.push_back(hist);
   }
   for (int i=0; i<nbins_centfwd; i++) {
      double centmin = bins_centfwd[i];
      double centmax = bins_centfwd[i+1];
      TH1F *hist = new TH1F(Form("hfwd_cent%.1f-%.1f",centmin,centmax),"",70,-0.15,0.2);
      hists_fwdcent.push_back(hist);
   }

   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      // if (ientry>1000) break;

      // skip the event if Gen_QQ_size != 1 for now
      TLorentzVector *tlvgenpl=NULL, *tlvgenmi=NULL, *tlvgenqq=NULL;
      if (!isdata) {
         b_Gen_QQ_size->GetEntry(ientry); //read only this branch
         if (Gen_QQ_size !=1) continue;

         // apply gen acceptance cuts
         b_Gen_QQ_mupl_4mom->GetEntry(ientry); //read only this branch
         b_Gen_QQ_mumi_4mom->GetEntry(ientry); //read only this branch
         tlvgenpl = (TLorentzVector*) Gen_QQ_mupl_4mom->At(0);
         tlvgenmi = (TLorentzVector*) Gen_QQ_mumi_4mom->At(0);
         if (!isGlobalMuonInAccept2015(tlvgenpl) || !isGlobalMuonInAccept2015(tlvgenmi)) continue;

         // fill denominators
         b_Gen_QQ_4mom->GetEntry(ientry); //read only this branch
         TLorentzVector *tlvgenqq = (TLorentzVector*) Gen_QQ_4mom->At(0);
         double genpt = tlvgenqq->Pt();
         b_Centrality->GetEntry(ientry); //read only this branch

         // check that the dimuon is inside the analysis bins
         bool gen_inbin = false;
         if (fabs(tlvgenqq->Rapidity())<1.6 && tlvgenqq->Pt()>6.5 && tlvgenqq->Pt()<30.) gen_inbin = true;
         if (fabs(tlvgenqq->Rapidity())>1.6 && fabs(tlvgenqq->Rapidity())<2.4 && tlvgenqq->Pt()>3. && tlvgenqq->Pt()<30.) gen_inbin = true;
         if (!gen_inbin) continue;
      }

      double weight = ispbpb ? fChain->GetWeight()*findNcoll(Centrality) : 1.;

      // is there at least one reco dimuon?
      b_Reco_QQ_size->GetEntry(ientry);
//.........这里部分代码省略.........
开发者ID:CMS-HIN-dilepton,项目名称:DimuonCADIs,代码行数:101,代码来源:ljpsiresol.C


注:本文中的TLorentzVector::Rapidity方法示例由纯净天空整理自Github/MSDocs等开源代码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。