本文整理汇总了C++中TLorentzVector::Vect方法的典型用法代码示例。如果您正苦于以下问题:C++ TLorentzVector::Vect方法的具体用法?C++ TLorentzVector::Vect怎么用?C++ TLorentzVector::Vect使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类TLorentzVector
的用法示例。
在下文中一共展示了TLorentzVector::Vect方法的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: CalcDeltaPhiRFRAME
// This is the pt corrected delta phi between the 2 leptons
// P and L2 are the 4-vectors for the 2 hemispheres, or in you case,
// the two leptons - setting mass to 0 should be fine
// MET is the MET 3 vector (don't forget to set the z-component of
// MET to 0)
// This function will do the correct Lorentz transformations of the
// leptons for you
double HWWKinematics::CalcDeltaPhiRFRAME(){
// first calculate pt-corrected MR
float mymrnew = CalcMRNEW();
// Now, boost lepton system to rest in z
// (approximate accounting for longitudinal boost)
TVector3 BL = L1.Vect()+L2.Vect();
BL.SetX(0.0);
BL.SetY(0.0);
BL = (1./(L1.P()+L2.P()))*BL;
L1.Boost(-BL);
L2.Boost(-BL);
// Next, calculate the transverse Lorentz transformation
// to go to Higgs approximate rest frame
TVector3 B = L1.Vect()+L2.Vect()+MET;
B.SetZ(0.0);
B = (-1./(sqrt(4.*mymrnew*mymrnew+B.Dot(B))))*B;
L1.Boost(B);
L2.Boost(B);
//Now, re-calculate the delta phi
// in the new reference frame:
return L1.DeltaPhi(L2);
}
示例2: CosThetaStar
double CosThetaStar(TLorentzVector p1, TLorentzVector p2){
TLorentzVector p = p1 + p2;
TVector3 theBoost = p.BoostVector();
TVector3 bostDir;
if ( theBoost.Mag() != 0 ) bostDir = theBoost.Unit(); // / theBoost.Mag());
else return -1;
p1.Boost(-theBoost);
if (p1.Vect().Mag()!=0) return p1.Vect().Dot(bostDir) / p1.Vect().Mag();
else return -1;
}
示例3: CalcMRNEW
// This is the pt corrected MR - here, we assume the pt
// of the Higgs is (MET+Pt+Qt);
// L1 and L2 are the 4-vectors for the 2 hemispheres, or in you case,
// the two leptons - setting mass to 0 should be fine
// MET is the MET 3 vector (don't forget to set the z-component of
// MET to 0)
// Also, 2 times this variable should give you the Higgs mass
double HWWKinematics::CalcMRNEW(){
TVector3 vI = MET+L1.Vect()+L2.Vect();
vI.SetZ(0.0);
double L1pL2 = CalcMR(); //Note - this calls the old MR function
double vptx = (L1+L2).Px();
double vpty = (L1+L2).Py();
TVector3 vpt;
vpt.SetXYZ(vptx,vpty,0.0);
float MR2 = 0.5*(L1pL2*L1pL2-vpt.Dot(vI)+L1pL2*sqrt(L1pL2*L1pL2+vI.Dot(vI)-2.*vI.Dot(vpt)));
return sqrt(MR2);
}
示例4:
Path::Path(const TLorentzVector& p4, const TVector3& origin, double field)
: m_unitDirection(p4.Vect().Unit()),
m_speed(p4.Beta() * gconstc),
m_origin(origin.X(), origin.Y(), origin.Z()),
m_field(field) {
m_points[papas::Position::kVertex] = m_origin;
}
示例5: CalcDeltaPhiNEW
// This is the pt corrected delta phi between the 2 mega-jets
// P and Q are the 4-vectors for the 2 hemispheres
// M is the MET 3 vector (don't forget to set the z-component of
// MET to 0)
// This function will do the correct Lorentz transformations of the
// leptons for you
double HWWKinematics::CalcDeltaPhiNEW(TLorentzVector P, TLorentzVector Q, TVector3 M){
// first calculate pt-corrected MR
float mymrnew = CalcMRNEW(L1,L2,MET);
//Next, calculate the transverse Lorentz transformation
TVector3 B = P.Vect()+Q.Vect()+MET;
B.SetZ(0.0);
B = (-1./(sqrt(4.*mymrnew*mymrnew+B.Dot(B))))*B;
P.Boost(B);
Q.Boost(B);
//Now, re-calculate the delta phi
// in the new reference frame:
return P.DeltaPhi(Q);
}
示例6: CalcUnboostedMTR
double HWWKinematics::CalcUnboostedMTR(TLorentzVector P, TLorentzVector Q, TVector3 M){
// first calculate pt-corrected MR
float mymrnew = CalcMRNEW(L1,L2,MET);
//Next, calculate the transverse Lorentz transformation
TVector3 B = P.Vect()+Q.Vect()+MET;
B.SetZ(0.0);
B = (-1./(sqrt(4.*mymrnew*mymrnew+B.Dot(B))))*B;
P.Boost(B);
Q.Boost(B);
//Now, re-calculate MTR in the new reference frame:
float mymtrnew = CalcMTRNEW(P, Q);
//R is now just the ratio of mymrnew and mymtrnew;
return mymtrnew;
}
示例7: cosThetaBoost
inline float kinematics::cosThetaBoost( TLorentzVector* pa, float ca,
TLorentzVector* pb, float cb )
{
// http://xrootd.slac.stanford.edu/BFROOT/www/doc/workbook_backup_010108/analysis/analysis.html
// A useful quantity in many analyses is the helicity angle.
// In the reaction Y -> X -> a + b, the helicity angle of
// particle a is the angle measured in the rest frame of the
//decaying parent particle, X, between the direction of the
// decay daughter a and the direction of the grandparent particle Y.
TLorentzVector pTmp = (*pa)+(*pb); // this is the mumu system (Z) 4vector
TVector3 ZboostVector = pTmp.BoostVector(); // this is the 3vector of the Z
TLorentzVector p; // this is the muon 4vector
if(ca<0) p.SetPxPyPzE(pa->Px(),pa->Py(),pa->Pz(),pa->E());
else if(cb<0) p.SetPxPyPzE(pb->Px(),pb->Py(),pb->Pz(),pb->E());
p.Boost( -ZboostVector ); // boost p to the dimuon CM (rest) frame
float cosThetaB = p.Vect()*pTmp.Vect()/(p.P()*pTmp.P());
//if (ySystem(pa,pb) < 0) cosThetaB *= -1.; // reclassify ???
return cosThetaB;
}
示例8: mu
////////////////////////////////////////////////////////////////////////
// Calculates theta and phi in HX frame
////////////////////////////////////////////////////////////////////////
pair<double, double> GetAngles_HX( TLorentzVector a, TLorentzVector b) {
TLorentzVector c = a+b; // JPsi momentum in lab frame
TVector3 bv = c.BoostVector();
TLorentzVector p1(0., 0., 1380., 1380.); // beam momentum in lab frame
TLorentzVector p2(0., 0., -1380., 1380.); // beam momentum in lab frame
p1.Boost(-bv);
p2.Boost(-bv);
TVector3 beam1 = p1.Vect().Unit(); // beam direction in JPsi rest frame
TVector3 beam2 = p2.Vect().Unit(); // beam direction in JPsi rest frame
TVector3 Z = c.Vect().Unit(); // JPsi direction in lab frame
TVector3 Y = beam1.Cross( beam2 ).Unit(); // the production plane normal
TVector3 X = Y.Cross(Z).Unit(); // completes the right-handed coordinate
a.Boost(-bv); // muon+ momentum in JPsi rest frame
TVector3 mu(a.Vect().Dot(X), a.Vect().Dot(Y), a.Vect().Dot(Z)); // transform to new coordinate
pair<double, double> angles;
angles.first = mu.Theta();
angles.second = mu.Phi()>0. ? mu.Phi() : mu.Phi()+2.*TMath::Pi();
return angles;
}
示例9: CalcDoubleDphiRFRAME
// This is the deltaphi between the di-lepton system and the Higgs
// boost in the approximate Higgs rest frame, or R-FRAME
// L1 and L2 are the 4-vectors for the 2 hemispheres, or in you case,
// the two leptons - setting mass to 0 should be fine
// MET is the MET 3 vector (don't forget to set the z-component of
// MET to 0)
// This function will do the correct Lorentz transformations of the
// leptons for you
double HWWKinematics::CalcDoubleDphiRFRAME(){
// first calculate pt-corrected MR
float mymrnew = CalcMRNEW();
TVector3 BL = L1.Vect()+L2.Vect();
BL.SetX(0.0);
BL.SetY(0.0);
BL = (1./(L1.P()+L2.P()))*BL;
L1.Boost(-BL);
L2.Boost(-BL);
//Next, calculate the transverse Lorentz transformation
TVector3 B = L1.Vect()+L2.Vect()+MET;
B.SetZ(0.0);
B = (-1./(sqrt(4.*mymrnew*mymrnew+B.Dot(B))))*B;
L1.Boost(B);
L2.Boost(B);
// Now, calculate the delta phi
// between di-lepton axis and boost
// in new reference frame
return B.DeltaPhi(L1.Vect()+L2.Vect());
}
示例10: Boost_To_Stop_Rest_Frame
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());
}
示例11: gravitonPythia
void gravitonPythia(){
gStyle->SetOptStat(0);
TreeReader data("pythia8_RS_WW.root");
TH1F* h_cosTh = new TH1F("h_cosTh","RS Graviton by Pythia8",100,-1,1);
h_cosTh->SetMinimum(0);
h_cosTh->Sumw2();
TH1F *h_cosThStar = (TH1F*)h_cosTh->Clone("h_cosThStar");
for(Long64_t ev = 0 ; ev < data.GetEntriesFast(); ev++){
data.GetEntry(ev);
Int_t nGenPar = data.GetInt("nGenPar");
Int_t* genParId = data.GetPtrInt("genParId");
Int_t* genParSt = data.GetPtrInt("genParSt");
Float_t* genParPt = data.GetPtrFloat("genParPt");
Float_t* genParEta = data.GetPtrFloat("genParEta");
Float_t* genParPhi = data.GetPtrFloat("genParPhi");
Float_t* genParM = data.GetPtrFloat("genParM");
// cos#theta_1 in the W rest frame
Int_t chgLepID = -1;
Int_t neuLepID = -1;
TLorentzVector chgLep(0,0,0,0);
TLorentzVector neuLep(0,0,0,0);
for(Int_t i = 0; i < nGenPar; i++){
if( genParSt[i] != 23 ) continue;
if( abs(genParId[i]) == 11 || abs(genParId[i]) == 13 || abs(genParId[i]) == 15 ) chgLepID = i;
if( abs(genParId[i]) == 12 || abs(genParId[i]) == 14 || abs(genParId[i]) == 16 ) neuLepID = i;
}
chgLep.SetPtEtaPhiM(genParPt[chgLepID],
genParEta[chgLepID],
genParPhi[chgLepID],
genParM[chgLepID]);
neuLep.SetPtEtaPhiM(genParPt[neuLepID],
genParEta[neuLepID],
genParPhi[neuLepID],
genParM[neuLepID]);
TLorentzVector Wb = chgLep + neuLep;
TVector3 WbP = Wb.Vect();
TVector3 bv = -Wb.BoostVector();
chgLep.Boost(bv);
TVector3 chgLepP = chgLep.Vect();
Double_t cosTh = TMath::Cos(chgLepP.Angle(WbP));
h_cosTh->Fill(cosTh);
// cos#theta* in the RSG rest frame
Int_t WplusID = -1;
Int_t WminusID = -1;
TLorentzVector Wplus(0,0,0,0);
TLorentzVector Wminus(0,0,0,0);
for(Int_t i = 0; i < nGenPar; i++){
if( genParSt[i] != 22 ) continue;
if( genParId[i] == +24 ) WplusID = i;
if( genParId[i] == -24 ) WminusID = i;
}
Wplus.SetPtEtaPhiM(genParPt[WplusID],
genParEta[WplusID],
genParPhi[WplusID],
genParM[WplusID]);
Wplus.SetPtEtaPhiM(genParPt[WminusID],
genParEta[WminusID],
genParPhi[WminusID],
genParM[WminusID]);
TLorentzVector RSG = Wplus + Wminus;
TVector3 RSGP = RSG.Vect();
TVector3 gv = -RSG.BoostVector();
Wplus.Boost(gv);
TVector3 WplusP = Wplus.Vect();
//.........这里部分代码省略.........
示例12: acceptance
void acceptance(const char* file = "upsilonGun.root", const char* outputName = "acceptance.root") {
Int_t genUpsSize;
Float_t genUpsPt[NMAX];
Float_t genUpsEta[NMAX];
Float_t genUpsPhi[NMAX];
Float_t genUpsRapidity[NMAX];
Int_t genMuSize;
Float_t genMuPt[NMAX];
Float_t genMuEta[NMAX];
Float_t genMuPhi[NMAX];
Int_t genMuCharge[NMAX];
Int_t recoMuSize;
Float_t recoMuPt[NMAX];
Float_t recoMuEta[NMAX];
Float_t recoMuPhi[NMAX];
Int_t recoMuCharge[NMAX];
TFile* f1 = new TFile(file);
TTree* t = (TTree*)f1->Get("UpsTree");
t->SetBranchAddress("genUpsSize",&genUpsSize);
t->SetBranchAddress("genUpsPt",genUpsPt);
t->SetBranchAddress("genUpsEta",genUpsEta);
t->SetBranchAddress("genUpsPhi",genUpsPhi);
t->SetBranchAddress("genUpsRapidity",genUpsRapidity);
t->SetBranchAddress("genMuSize",&genMuSize);
t->SetBranchAddress("genMuPt",genMuPt);
t->SetBranchAddress("genMuEta",genMuEta);
t->SetBranchAddress("genMuPhi",genMuPhi);
t->SetBranchAddress("genMuCharge",genMuCharge);
t->SetBranchAddress("recoMuSize",&recoMuSize);
t->SetBranchAddress("recoMuPt",recoMuPt);
t->SetBranchAddress("recoMuEta",recoMuEta);
t->SetBranchAddress("recoMuPhi",recoMuPhi);
t->SetBranchAddress("recoMuCharge",recoMuCharge);
TH2F* genPtRap = new TH2F("genUps","",4,0,20,4,-2,2);
TH1F* genPt = new TH1F("genUps1D","",4,0,20);
genPtRap->Sumw2();
genPt->Sumw2();
genPtRap->SetTitle(";Upsilon pT (GeV/c);Upsilon Rapidity;");
genPt->SetTitle(";Upsilon pT (GeV/c);Acceptance;");
TH2F* recoPtRap = (TH2F*)genPtRap->Clone("recoUps");
TH1F* recoPt = (TH1F*)genPt->Clone("recoUps1D");
recoPtRap->Sumw2();
recoPt->Sumw2();
for(int i=0; i<t->GetEntries(); i++){
if(i%100000 == 0)
std::cout<<i<<std::endl;
t->GetEntry(i);
// calculate cosTheta
TLorentzVector genUps; genUps.SetPtEtaPhiM(genUpsPt[0], genUpsEta[0], genUpsPhi[0], 9.46);
TLorentzRotation boost(-genUps.BoostVector());
int mp = genMuCharge[0]>0 ? 0 : 1;
TLorentzVector genMuPlus; genMuPlus.SetPtEtaPhiM(genMuPt[mp], genMuEta[mp], genMuPhi[mp], 0.106);
genMuPlus *= boost;
Float_t cosThetaStar = genMuPlus.Vect().Dot(genUps.Vect())/genMuPlus.Vect().Mag()/genUps.Vect().Mag();
// set the weight
Float_t weight = 1;
genPtRap->Fill( genUpsPt[0], genUpsRapidity[0], weight );
genPt->Fill( genUpsPt[0], weight );
Float_t recoUpsPt = 0;
Float_t recoUpsRapidity = 0;
double minDeltaM = 1000;
for(int tr1=0; tr1<recoMuSize; tr1++){
for(int tr2=tr1+1; tr2<recoMuSize; tr2++){
if ( recoMuCharge[tr1]*recoMuCharge[tr2] == -1 && recoMuPt[tr1] >= 3.5 && recoMuPt[tr2] >= 3.5 && fabs(recoMuEta[tr1]) < 2.1 && fabs(recoMuEta[tr2]) < 2.1 ){
TLorentzVector mu1; mu1.SetPtEtaPhiM(recoMuPt[tr1], recoMuEta[tr1], recoMuPhi[tr1], 0.106);
TLorentzVector mu2; mu2.SetPtEtaPhiM(recoMuPt[tr2], recoMuEta[tr2], recoMuPhi[tr2], 0.106);
TLorentzVector recoUps(mu1 + mu2);
double deltaM = fabs(recoUps.M()-9.46);
if( deltaM < minDeltaM ){
recoUpsPt = recoUps.Pt();
recoUpsRapidity = recoUps.Rapidity();
minDeltaM = deltaM;
}
}
}
}
if( minDeltaM < 1.0 ){
recoPtRap->Fill( recoUpsPt, recoUpsRapidity, weight );
recoPt->Fill( recoUpsPt, weight );
}
}
TFile out(outputName,"recreate");
TH2F* acc = (TH2F*)genPtRap->Clone("acceptance");
TH1F* acc1D = (TH1F*)genPt->Clone("acceptance1D");
acc->Sumw2();
acc1D->Sumw2();
acc->Divide(recoPtRap,genPtRap,1,1,"B");
acc1D->Divide(recoPt,genPt,1,1,"B");
acc->Write();
acc1D->Write();
out.Close();
}
示例13: main
int main (int argc, char ** argv)
{
if(argc < 3)
{
cout << "Usage: " << argv[0]
<< " input.lhe output.lhe" << endl ;
return -1;
}
std::ifstream ifs (argv[1]) ;
LHEF::Reader reader (ifs) ;
ofstream outputStream (argv[2]) ;
LHEF::Writer writer (outputStream) ;
writer.headerBlock () << reader.headerBlock ;
writer.initComments () << reader.initComments ;
writer.heprup = reader.heprup ;
writer.init () ;
//PG mu massless in phantom
// float k2 = 0.1056583715 * 0.1056583715 - 1.77682 * 1.77682 ; // GeV -3.14592562093
float k2 = 0. - 1.77682 * 1.77682 ; // GeV -3.14592562093
int count = 0 ;
//PG loop over input events
while (reader.readEvent ())
{
++count ;
if ( reader.outsideBlock.length ()) std::cout << reader.outsideBlock;
// loop over particles in the event
for (int iPart = 0 ; iPart < reader.hepeup.IDUP.size (); ++iPart)
{
// outgoing particles
if (reader.hepeup.ISTUP.at (iPart) == 1)
{
if (abs (reader.hepeup.IDUP.at (iPart)) == 13)
{
TLorentzVector dummy
(
reader.hepeup.PUP.at (iPart).at (0), // px
reader.hepeup.PUP.at (iPart).at (1), // py
reader.hepeup.PUP.at (iPart).at (2), // pz
reader.hepeup.PUP.at (iPart).at (3) // E
) ;
float p2 = dummy.Vect ().Mag2 () ;
float scale = sqrt (1 + k2 / p2) ;
if (p2 < (-1 * k2))
{
cout << "warning: p2 is smaller than the mass difference " << p2 << endl ;
scale = 1 ;
}
reader.hepeup.PUP.at (iPart).at (0) *= scale ; // px
reader.hepeup.PUP.at (iPart).at (1) *= scale ; // px
reader.hepeup.PUP.at (iPart).at (2) *= scale ; // px
if (reader.hepeup.IDUP.at (iPart) == 13) reader.hepeup.IDUP.at (iPart) = 15 ;
if (reader.hepeup.IDUP.at (iPart) == -13) reader.hepeup.IDUP.at (iPart) = -15 ;
}
if (reader.hepeup.IDUP.at (iPart) == 14) reader.hepeup.IDUP.at (iPart) = 16 ;
if (reader.hepeup.IDUP.at (iPart) == -14) reader.hepeup.IDUP.at (iPart) = -16 ;
} // outgoing particles
} // loop over particles in the event
writer.eventComments () << reader.eventComments ;
writer.hepeup = reader.hepeup ;
bool written = writer.writeEvent () ;
if (!written)
{
cout << "warning: event " << count << " not written" << endl ;
}
} //PG loop over input events
cout << "end loop over " << count << " events" << endl ;
return 0 ;
}
示例14: uDstChain
void
fillUdstDataIntoMassBins_example(const string& inFileNamePattern = "fillUdstDataIntoMassBins_example.root",
const string& dirName = "./test",
const long int maxNmbEvents = -1,
const unsigned int nmbMassBins = 50,
const double massBinWidth = 40, // [MeV/c^2]
const double massRangeMin = 500, // [MeV/c^2]
const string& uDstTreeName = "pwaDataTree",
const string& pwaTreeName = "rootPwaEvtTree",
const long int treeCacheSize = 25000000, // 25 MByte ROOT tree read cache
const bool debug = false)
{
const string prodKinPartNamesObjName = "prodKinParticles";
const string prodKinMomentaLeafName = "prodKinMomenta";
const string decayKinPartNamesObjName = "decayKinParticles";
const string decayKinMomentaLeafName = "decayKinMomenta";
TStopwatch timer;
timer.Start();
printInfo << "reading uDST file(s) '" << inFileNamePattern << "'" << endl
<< " writing " << nmbMassBins << " mass bins in mass interval "
<< "[" << massRangeMin << ", " << massRangeMin + nmbMassBins * massBinWidth << "] "
<< "MeV/c^2 to '" << dirName << "'" << endl
<< " reading uDST data from tree '" << uDstTreeName << "'" << endl
<< " writing PWA data to tree '" << pwaTreeName << "'" << endl;
// create chain and connect tree leaf variables to branches
TChain uDstChain(uDstTreeName.c_str());
if (uDstChain.Add(inFileNamePattern.c_str()) < 1)
printWarn << "no events in uDST input file(s) '" << inFileNamePattern << "'" << endl;
const long int nmbEventsUdstChain = uDstChain.GetEntries();
uDstChain.GetListOfFiles()->ls();
// !!! <channel-dependent part> !!!
// connect tree leafs
TLorentzVector* photons[4] = {0, 0, 0, 0};
TLorentzVector* piMinus = 0;
TLorentzVector* beam = 0;
TLorentzVector* recoil = 0;
uDstChain.SetBranchAddress("gamma1", &(photons[0]));
uDstChain.SetBranchAddress("gamma2", &(photons[1]));
uDstChain.SetBranchAddress("gamma3", &(photons[2]));
uDstChain.SetBranchAddress("gamma4", &(photons[3]));
uDstChain.SetBranchAddress("pi_out", &piMinus);
uDstChain.SetBranchAddress("pi_in", &beam);
uDstChain.SetBranchAddress("proton", &recoil);
uDstChain.SetCacheSize(treeCacheSize);
uDstChain.AddBranchToCache("gamma1", true);
uDstChain.AddBranchToCache("gamma2", true);
uDstChain.AddBranchToCache("gamma3", true);
uDstChain.AddBranchToCache("gamma4", true);
uDstChain.AddBranchToCache("pi_out", true);
uDstChain.AddBranchToCache("pi_in", true);
uDstChain.AddBranchToCache("proton", true);
uDstChain.StopCacheLearningPhase();
// !!! </channel-dependent part> !!!
// create directories and .root files
vector<TFile*> pwaDataFiles;
vector<TTree*> pwaDataTrees;
if (not createMassBinFiles(pwaDataFiles, pwaDataTrees, dirName, nmbMassBins, massBinWidth,
massRangeMin, pwaTreeName)) {
printErr << "there were problems creating the mass bin directories/files. aborting." << endl;
return;
}
printSucc << "created " << pwaDataFiles.size() << " directories/files" << endl;
// write arrays with production and decay particle names to root files
{
TClonesArray prodKinPartNames ("TObjString", 1);
TClonesArray decayKinPartNames("TObjString", 3);
// !!! <channel-dependent part> !!!
new (prodKinPartNames [0]) TObjString("pi-"); // beam particle
new (decayKinPartNames[0]) TObjString("pi0");
new (decayKinPartNames[1]) TObjString("pi0");
new (decayKinPartNames[2]) TObjString("pi-");
// !!! </channel-dependent part> !!!
for (unsigned int i = 0; i < pwaDataFiles.size(); ++i) {
pwaDataFiles[i]->cd();
prodKinPartNames.Write (prodKinPartNamesObjName.c_str (), TObject::kSingleKey);
decayKinPartNames.Write(decayKinPartNamesObjName.c_str(), TObject::kSingleKey);
}
printSucc << "wrote particle name arrays to all files. "
<< "beam = 'pi-', decay = {'pi0', 'pi0', 'pi-'}." << endl;
}
// create tree leafs
{
TClonesArray* prodKinMomenta = new TClonesArray("TVector3");
TClonesArray* decayKinMomenta = new TClonesArray("TVector3");
const int splitLevel = 99;
const int bufSize = 256000;
for (unsigned int i = 0; i < pwaDataTrees.size(); ++i) {
pwaDataTrees[i]->Branch(prodKinMomentaLeafName.c_str(), "TClonesArray", &prodKinMomenta,
bufSize, splitLevel);
pwaDataTrees[i]->Branch(decayKinMomentaLeafName.c_str(), "TClonesArray", &decayKinMomenta,
bufSize, splitLevel);
//.........这里部分代码省略.........
示例15: TClonesArray
// fills event data into correct mass bin
bool
writeEvent(vector<TTree*>& pwaTrees,
const TLorentzVector& beamLv,
// !!! <channel-dependent part> !!!
const TLorentzVector& piZero0,
const TLorentzVector& piZero1,
const TLorentzVector& piMinus,
// !!! </channel-dependent part> !!!
const double XMass, // [GeV/c^2]
const unsigned int nmbMassBins = 50,
const double massBinWidth = 50, // [MeV/c^2]
const double massRangeMin = 500, // [MeV/c^2]
const string& prodKinMomentaLeafName = "prodKinMomenta",
const string& decayKinMomentaLeafName = "decayKinMomenta",
const bool debug = false)
{
const double mass = 1000 * XMass; // convert from GeV/c^2 to MeV/c^2
// make sure that mass is in range
if ((mass < massRangeMin) or (mass > (massRangeMin + nmbMassBins * massBinWidth)))
return false;
const unsigned int bin = (unsigned int) ((mass - massRangeMin) / massBinWidth);
if (not pwaTrees[bin]) {
printWarn << "null pointer for tree for mass bin [" << massRangeMin + bin * massBinWidth << ", "
<< massRangeMin + (bin + 1) * massBinWidth << "]" << endl;
return false;
}
// fill tree
if (debug)
printDebug << "filling tree for bin " << bin << " = ["
<< massRangeMin + bin * massBinWidth << ", "
<< massRangeMin + (bin + 1) * massBinWidth << "] MeV/c^2" << endl;
// create tree leafs
static TClonesArray* prodKinMomenta = new TClonesArray("TVector3");
static TClonesArray* decayKinMomenta = new TClonesArray("TVector3");
// connect leaf variables to tree branches or create branches, if they don't exist yet
TTree* outTree = pwaTrees[bin];
if (outTree->SetBranchAddress(prodKinMomentaLeafName.c_str(), &prodKinMomenta) < 0) {
printWarn << "could not connect variable to branch '" << prodKinMomentaLeafName << "'. "
<< "skipping." << endl;
return false;
}
if (outTree->SetBranchAddress(decayKinMomentaLeafName.c_str(), &decayKinMomenta) < 0) {
printWarn << "could not connect variable to branch '" << prodKinMomentaLeafName << "'. "
<< "skipping." << endl;
return false;
}
// clear arrays
prodKinMomenta->Clear ();
decayKinMomenta->Clear();
// set leaf variables
// beam particle
new ((*prodKinMomenta)[0]) TVector3(beamLv.Vect());
// !!! <channel-dependent part> !!!
// for target particle elastic scattering is assumed
// outgoing hadrons
new ((*decayKinMomenta)[0]) TVector3(piZero0.Vect());
new ((*decayKinMomenta)[1]) TVector3(piZero1.Vect());
new ((*decayKinMomenta)[2]) TVector3(piMinus.Vect());
// !!! </channel-dependent part> !!!
// fill tree
outTree->Fill();
return true;
}