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

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;
}

void rochcor2012::musclefit_data( TLorentzVector& mu, TLorentzVector& mubar){

  float dpar1 = 0.0;
  float dpar2 = 0.0;
  float epar1 = 0.0;
  float epar2 = 0.0;
  
  if(fabs(mu.PseudoRapidity())<=0.9){
    dpar1 = d0par;
    epar1 = e0par;
  }else if(mu.PseudoRapidity()>0.9){
    dpar1 = d1par;
    epar1 = e1par;
  }else if(mu.PseudoRapidity()<-0.9){
    dpar1 = d2par;
    epar1 = e2par;
  }

  if(fabs(mubar.PseudoRapidity())<=0.9){
    dpar2 = d0par;
    epar2 = e0par;
  }else if(mubar.PseudoRapidity()>0.9){
    dpar2 = d1par;
    epar2 = e1par;
  }else if(mubar.PseudoRapidity()<-0.9){
    dpar2 = d2par;
    epar2 = e2par;
  }

  float corr1 = 1.0 + bpar*mu.Pt() + (-1.0)*cpar*mu.Pt()*TMath::Sign(float(1.0),float(mu.PseudoRapidity()))*TMath::Power(mu.PseudoRapidity(),2)
    + (-1.0)*dpar1*mu.Pt()*sin(mu.Phi() + epar1);
  float corr2 = 1.0 + bpar*mubar.Pt() + (1.0)*cpar*mubar.Pt()*TMath::Sign(float(1.0),float(mubar.PseudoRapidity()))*TMath::Power(mubar.PseudoRapidity(),2)
    + (1.0)*dpar2*mubar.Pt()*sin(mubar.Phi() + epar2);
  
  float px1 = mu.Px();
  float py1 = mu.Py();
  float pz1 = mu.Pz();
  float e1 = mu.E();
  
  float px2 = mubar.Px();
  float py2 = mubar.Py();
  float pz2 = mubar.Pz();
  float e2 = mubar.E();

  px1 *= corr1;
  py1 *= corr1;
  pz1 *= corr1;
  e1 *= corr1;
  
  px2 *= corr2;
  py2 *= corr2;
  pz2 *= corr2;
  e2 *= corr2;
  
  mu.SetPxPyPzE(px1,py1,pz1,e1);
  mubar.SetPxPyPzE(px2,py2,pz2,e2);

}

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;

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

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);
}

//==============================================
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;
//.........这里部分代码省略.........