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

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

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

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

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

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

}

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

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

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

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

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

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

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

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

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

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

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