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

double Rho0Omega::Getp( const double& mpipi ) const {
  const TLorentzVector p4pipi =
    (BW_BW::daughterP4(0)+BW_BW::daughterP4(1))*(1.0/GeV);
  const TLorentzVector p4init =
    _theDecay.getParent()->getVal().p(getEvent())*(1.0/GeV);
  const TLorentzVector p4rec = p4init-p4pipi;

  const double minit = p4init.M();
  const double mrec = p4rec.M();

  const double m2init = minit*minit;
  const double mpipi_p2_mrec = (mpipi + mrec)*(mpipi + mrec);
  const double mpipi_m2_mrec = (mpipi - mrec)*(mpipi - mrec);

  const double p =
    sqrt((m2init-mpipi_p2_mrec)*(m2init-mpipi_m2_mrec))/
    (2.0*minit);

  if( (m2init-mpipi_p2_mrec)*(m2init-mpipi_m2_mrec) < 0.0 ){
    std::cout << "p is negative. This should never happen." << std::endl;
    exit(1);
  }

  return p;
}

// Takes as input a given process with corresponding luminosity, as well as sample CME, cuts, directory location, and output file name - saves corresponding histogram in root file.
void processtree(TString sample, Double_t weight, TLorentzVector CME, TString filename, TString directory, TString param, TString cut, std::clock_t start, Int_t all){

    //Open File
    TFile *f = new TFile(directory);
    TTree *ttree = (TTree*)f->Get(sample);

    TString recohistname, paramhistname;
    if (all==0){
        recohistname = "_recoilmassHIST";
        paramhistname = "_parameterHIST";
    }
    if (all==1){
        recohistname = "_recoilmassHIST_all";
        paramhistname = "_parameterHIST_all";
    }

    TString histname = sample + recohistname;
    TH1D *recoilmassHIST   = new TH1D(histname,histname,CME.M()/2,0,CME.M());
    //recoilmassHIST->SetCanExtend(kAllAxes);
    ttree->Project(histname, "recoilmass", cut);

    histname = sample + paramhistname;
    TH1D *parameterHIST   = new TH1D(histname,histname,CME.M()/2,1,-1);
    //parameterHIST->SetCanExtend(kAllAxes);
    ttree->Project(histname, param, cut);

    //Normalize Hist
    recoilmassHIST->Scale(1/recoilmassHIST->GetMaximum());
    parameterHIST->Scale(1/parameterHIST->GetMaximum()) ;

    //Normalize Hist
//    recoilmassHIST->Scale(weight);
//    parameterHIST->Scale(weight) ;

    // Writes yield of each histogram for each cut to text file
    ofstream myfile;
    Double_t yield = recoilmassHIST->Integral(); 
    myfile.open ("yield.txt", ios::app);
        cout << "Yield for " << sample << " for Cut " << param << " is " << yield << endl;
        myfile << param << " " << all << " " << sample << " " << yield << "\n"; //write to file
    myfile.close();

    if (recoilmassHIST->Integral() == 0){
        recoilmassHIST->Fill(0);
    }
    if (parameterHIST->Integral() == 0){
        parameterHIST->Fill(0);
    }

    TFile g(filename, "update");
    recoilmassHIST->Write();
    parameterHIST->Write();
    g.Close();
}

void fillCorrectMassBin(TH1D* numhist[nMassBins], TLorentzVector sumVect, double phiSR, TH1D* polHist[nMassBins][nPhisrBins], double polarization, TH1D* massHist[nMassBins]){
	
	for (int iMass = 0; iMass < nMassBins; iMass++)
	{
		if (sumVect.M() >= massBinStart[iMass] && sumVect.M() < massBinEnd[iMass])
		{
			numhist[iMass]->Fill(phiSR);
			massHist[iMass]->Fill(sumVect.M());

			int phiSRbin = numhist[iMass]->FindBin(phiSR);
			polHist[iMass][phiSRbin-1]->Fill(polarization);
		}
	}
}

void smearEmEnergy( const EnergyScaleCorrection_class& egcor, TLorentzVector& p, const ZGTree& myTree, float r9 ) {

    float smearSigma =  egcor.getSmearingSigma(myTree.run, fabs(p.Eta())<1.479, r9, p.Eta(), p.Pt(), 0., 0.);
    float cor = myRandom_.Gaus( 1., smearSigma );
    p.SetPtEtaPhiM( p.Pt()*cor, p.Eta(), p.Phi(), p.M() );

}

void ToyEvent7::GenerateSignalKinematics(TRandom *rnd,Bool_t isData) {

   // fake decay of Z0 to two fermions
   double M0=91.1876;
   double Gamma=2.4952;
   // generated mass
   do {
      fMGen[2]=rnd->BreitWigner(M0,Gamma);
   } while(fMGen[2]<=0.0);

   double N_ETA=3.0;
   double MU_PT=5.;
   double SIGMA_PT=2.0;
   double DECAY_A=0.2;
   if(isData) {
      //N_ETA=2.5;
      MU_PT=6.;
      SIGMA_PT=1.8;
      //DECAY_A=0.5;
   }
   fEtaGen[2]=TMath::Power(rnd->Uniform(0,1.5),N_ETA);
   if(rnd->Uniform(-1.,1.)<0.) fEtaGen[2] *= -1.;
   fPhiGen[2]=rnd->Uniform(-M_PI,M_PI);
   do {
      fPtGen[2]=rnd->Landau(MU_PT,SIGMA_PT);
   } while((fPtGen[2]<=0.0)||(fPtGen[2]>500.));
   //========================== decay
   TLorentzVector sum;
   sum.SetPtEtaPhiM(fPtGen[2],fEtaGen[2],fPhiGen[2],fMGen[2]);
   // boost into lab-frame
   TVector3 boost=sum.BoostVector();
   // decay in rest-frame

   TLorentzVector p[3];
   double m=MASS1;
   double costh;
   do {
      double r=rnd->Uniform(-1.,1.);
      costh=r*(1.+DECAY_A*r*r);
   } while(fabs(costh)>=1.0);
   double phi=rnd->Uniform(-M_PI,M_PI);
   double e=0.5*sum.M();
   double ptot=TMath::Sqrt(e+m)*TMath::Sqrt(e-m);
   double pz=ptot*costh;
   double pt=TMath::Sqrt(ptot+pz)*TMath::Sqrt(ptot-pz);
   double px=pt*cos(phi);
   double py=pt*sin(phi);
   p[0].SetXYZT(px,py,pz,e);
   p[1].SetXYZT(-px,-py,-pz,e);
   for(int i=0;i<2;i++) {
      p[i].Boost(boost);
   }
   p[2]=p[0]+p[1];
   for(int i=0;i<3;i++) {
      fPtGen[i]=p[i].Pt();
      fEtaGen[i]=p[i].Eta();
      fPhiGen[i]=p[i].Phi();
      fMGen[i]=p[i].M();
   }
}

TLorentzVector smearMom(TLorentzVector const& b,TF1 const * const fMomResolution)
{
  float const pt = b.Perp();
  float const sPt = gRandom->Gaus(pt,pt*fMomResolution->Eval(pt));

  TLorentzVector sMom;
  sMom.SetXYZM(sPt*cos(b.Phi()),sPt*sin(b.Phi()),sPt*sinh(b.PseudoRapidity()),b.M());
  return sMom;
}

double InvMass(double pt1, double eta1, double phi1, double m1, double pt2, double eta2, double phi2, double m2){

	TLorentzVector vec1;
	TLorentzVector vec2;
	vec1.SetPtEtaPhiM(pt1, eta1, phi1, m1);
	vec2.SetPtEtaPhiM(pt2, eta2, phi2, m2);
	TLorentzVector Z = vec1 + vec2;
	return  Z.M();

}

Double_t EventClassifier::ScoreEvent(struct WWScatteringComponents *event) {
    TLorentzVector *jetPair = new TLorentzVector(*event->positiveJet + *event->negativeJet);
    HadronicJetAbsEta = abs(event->hadronicJet->Eta());
    HadronicJetPT = event->hadronicJet->Phi();
    HadronicJetMass = event->hadronicJet->M();
    MissingET = event->missingET->E();
    Mjj = jetPair->M();
    LeptonAbsEta = abs(event->lepton->Eta());
    LeptonPT = event->lepton->Pt();

    return reader->EvaluateMVA("LD");
}

//these can be done better, pass in a specific struct and populate that struct
//could reduce these 6 calls into 1 or 2
void DalitzChiSq2::setap2(TLorentzVector vec){
	ap2.v.SetXYZM(vec.Px(),vec.Py(),vec.Pz(),vec.M());

	double Pt = TMath::Sqrt(vec.Px()*vec.Px() + vec.Py()*vec.Py());
	double p = TMath::Sqrt(Pt*Pt + vec.Pz()*vec.Pz());
	double thet = TMath::ACos(vec.Pz()/p);
	
	ap2.pt=Pt;
	ap2.theta=thet;
	ap2.x_m=1/Pt;
	
}

void Boost_To_Stop_Rest_Frame(TLorentzVector& stop4, TLorentzVector& chargino4, TLorentzVector& b4, TLorentzVector& neutralino4, TLorentzVector& W4, TLorentzVector& up4, TLorentzVector& down4, TLorentzVector& s4)
{
    TVector3 betaV(-stop4.Px()/stop4.Energy(),-stop4.Py()/stop4.Energy(),-stop4.Pz()/stop4.Energy());
    stop4.Boost(betaV);
    chargino4.Boost(betaV);
    b4.Boost(betaV);
    neutralino4.Boost(betaV);
    W4.Boost(betaV);
    up4.Boost(betaV);
    down4.Boost(betaV);
    s4.SetE(chargino4.P()/chargino4.M());
    s4.SetVect(chargino4.Vect().Unit()*chargino4.Gamma());
}

Double_t ACAna::Minv(TLorentzVector *particle[],Int_t Number) {
  TLorentzVector *sparticle = new TLorentzVector(0,0,0,0);
  for (Int_t i=0;i<Number;++i) { 
      //cout<<"Energie: "<< particle[i]->E()<<" i: "<<i<<endl;
    *sparticle += *(particle[i]); 
  }
  Double_t Minv =sparticle->M(); 
  delete sparticle;
  return Minv;
  //Double_t Minvsq = sE*sE - sPx*sPx - sPy*sPy- sPz*sPz;
  //return Minvsq >= 0. ? sqrt(Minvsq) : -sqrt(-Minvsq);

}

void testTDime(Int_t nev = 100) {

  gSystem->Load("libEVGEN");
  gSystem->Load("libTDime");
  gSystem->Load("libdime");

  TDime* dime = new TDime();
  dime->SetEnergyCMS(7000.0);
  dime->SetYRange(-2.0, 2.0);   // Set rapidity range of mesons
  dime->SetMinPt(0.1);          // Minimum pT of mesons
  dime->Initialize();

  // (pi+pi-) histograms
  TH1* hM = new TH1D("hM", "DIME #pi^{+}#pi^{-};M_{#pi^{+}#pi^{-}} #[]{GeV/#it{c}^{2}}", 100,  0.0, 5.0);

  TClonesArray* particles = new TClonesArray("TParticle", 6);
  TParticle* part = NULL;
  TLorentzVector v[2];
  TLorentzVector vSum;

  // Event loop
  for (Int_t i = 0; i < nev; ++i) {

    dime->GenerateEvent();
    Int_t np = dime->ImportParticles(particles, "All");
    printf("\n DIME Event %d: Imported %3d particles \n", i, np);

    Int_t nPrimary = 0;

    // Loop over pion (j = 4,5) tracks
    for (Int_t j = 4; j < 6; ++j) {
      part = (TParticle*) particles->At(j); // Choose the particle
      part->Print();
      part->Momentum(v[nPrimary]);          // Copy content to v
      nPrimary++;
    }
    //particles.Clear();

    // 4-vector sum
    vSum = v[0] + v[1];

    // Fill pi+pi- histograms
    hM->Fill(vSum.M());
  }

  // Save plots as pdf
  hM->Draw();    c1->SaveAs("massTDime.pdf");

}

float getMT(TLorentzVector pZ, TLorentzVector pH)
{
  //take MET from pZ
  float myMET = pZ.Pt();
  float myMETx = pZ.Px();
  float myMETy = pZ.Py();

  float Et = pH.Et();
  float px = pH.Px();
  float py = pH.Py();
  float m  = pH.M();

  float MT = -99;
  float MT2 = m*m + 2*( Et * myMET - (px*myMETx + py*myMETy) ); 
  if(MT2>0.) {MT=sqrt(MT2);}
  
  return MT;
}

void smearEmEnergy( TLorentzVector& p ) {

    float smearEBlowEta     = 0.013898;
    float smearEBhighEta    = 0.0189895;
    float smearEElowEta     = 0.027686;
    float smearEEhighEta    = 0.031312;
    float theGaussMean      = 1.;
    float pt = p.Pt();
    float eta = p.Eta();
    float phi = p.Phi();
    float mass = p.M();
    float theSmear = 0.;
    if      (fabs(eta)<1.                   )  theSmear = smearEBlowEta;
    else if (fabs(eta)>=1.  && fabs(eta)<1.5)  theSmear = smearEBhighEta;
    else if (fabs(eta)>=1.5 && fabs(eta)<2. )  theSmear = smearEElowEta;
    else if (fabs(eta)>=2.  && fabs(eta)<2.5)  theSmear = smearEEhighEta;
    float fromGauss = myRandom_.Gaus(theGaussMean,theSmear);
    p.SetPtEtaPhiM( fromGauss*pt, eta, phi, mass ); // keep mass and direction same

}

void ptLooper::doLoop() {

  for (int i = ievent_; i < max_entry; i++) {
    if (i % 100000 == 0) cout << i << " / " << max_entry << endl;
    tree->GetEntry(i);

    TLorentzVector muPN;
    muPN = *muonP_p4 + *muonN_p4;
    muPN_mass_h.Fill(muPN.M());

    mass_h.Fill(dimuon_p4->M());

    dimuonPt_h.Fill(dimuon_p4->Pt());
    dimuonEta_h.Fill(dimuon_p4->Eta());

    dimuonpteta_h.Fill(dimuon_p4->Eta(), dimuon_p4->Pt());

    muonpteta_h.Fill(muonP_p4->Eta(), muonP_p4->Pt());
    muonpteta_h.Fill(muonN_p4->Eta(), muonN_p4->Pt());

    if (muonP_p4->Pt() > 4. && muonN_p4->Pt()>4.) {
      dimuonPt4_h.Fill(dimuon_p4->Pt());
      dimuonEta4_h.Fill(dimuon_p4->Eta());
    }

    if (simul_Jpsi_trigger(muonP_p4, muonN_p4, dimuon_p4)) {
      dimuonPtTriggered_h.Fill(dimuon_p4->Pt());
    }
  }

  histos1D.push_back(dimuonPt_h);
  histos1D.push_back(dimuonPt4_h);
  histos1D.push_back(dimuonEta_h);
  histos1D.push_back(dimuonEta4_h);
  histos1D.push_back(muPN_mass_h);
  histos1D.push_back(mass_h);
  histos1D.push_back(dimuonPtTriggered_h);

}