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

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

inline TLorentzVector* kinematics::Boost( TLorentzVector* pa, TLorentzVector* pb, TLorentzVector* p )
{
	TLorentzVector pTmp = (*pa)+(*pb);
	TLorentzVector* pBoosted = (TLorentzVector*)p->Clone("");
	pBoosted->Boost(-1.*pTmp.BoostVector());
	return pBoosted;
}

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

void Tran_to_fixed_target (TLorentzVector e_vec, TLorentzVector d_vec,
                            TLorentzVector &e_ft, TLorentzVector  &d_ft){
   //input: 
   //   e_vec - electron 4 vector in some frame
   //   d_vec - ion 4 vector in some frame
   //output:
   //   e_ft - electron 4 vector in rest frame of ion (fixed target frame)
   //   d_ft - ion 4 vector in its rest frame 
   e_ft = e_vec;
   d_ft = d_vec;
   
   d_ft.Boost(-d_vec.BoostVector()); //should be (0,0,0, mass)
   e_ft.Boost(-d_vec.BoostVector());
   return;


}

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

void analysis_pbarp_Xi_test(int nevts=0){
  
  TDatabasePDG::Instance()-> AddParticle("pbarpSystem","pbarpSystem", 1.9, kFALSE, 0.1, 0,"", 88888);
  
  TStopwatch timer; 



  //Output File
  TString Path = "/private/puetz/mysimulations/analysis/pbarp_Xiplus_Ximinus/idealtracking/10000_events/";
  TString outPath = Path;
  TString OutputFile = outPath + "analysis_output_test.root";
  
  //Input simulation Files
  TString inPIDFile = Path + "pid_complete.root";
  TString inParFile = Path + "simparams.root";
  TString PIDParFile = TString( gSystem->Getenv("VMCWORKDIR")) + "/macro/params/all.par";
  
  //Initialization
  FairLogger::GetLogger()->SetLogToFile(kFALSE);
  FairRunAna* RunAna = new FairRunAna();
  FairRuntimeDb* rtdb = RunAna->GetRuntimeDb();
  RunAna->SetInputFile(inPIDFile);
  
  //setup parameter database
  FairParRootFileIo* parIo = new FairParRootFileIo();
  parIo->open(inParFile);
  FairParAsciiFileIo* parIoPID = new FairParAsciiFileIo();
  parIoPID->open(PIDParFile.Data(),"in");
  
  rtdb->setFirstInput(parIo);
  rtdb->setSecondInput(parIoPID);
  rtdb->setOutput(parIo);

  RunAna->SetOutputFile(OutputFile);
  RunAna->Init();

  //*** create tuples
  RhoTuple * ntpMC = new RhoTuple("ntpMC", "MCTruth info");
  RhoTuple * ntpPiMinus = new RhoTuple("ntpPiMinus", "PiMinus info");
  RhoTuple * ntpPiPlus = new RhoTuple("ntpPiPlus", "PiPlus info");
  RhoTuple * ntpProton = new RhoTuple("ntpProton", "Proton info");
  RhoTuple * ntpAntiProton = new RhoTuple("ntpAntiProton", "Antiproton info");
  RhoTuple * ntpLambda0 = new RhoTuple("ntpLambda0", "Lambda0 info");
  RhoTuple * ntpAntiLambda0 = new RhoTuple("ntpAntiLambda0", "AntiLambda0 info");
  RhoTuple * ntpXiMinus = new RhoTuple("ntpXiMinus", "XiMinus info");
  RhoTuple * ntpXiPlus = new RhoTuple("ntpXiPlus", "XiPlus info");
  RhoTuple * ntpXiSys = new RhoTuple("ntpXiSys", "XiMinus XiPlus system info");

  //Create output file 
  TFile *out = TFile::Open(outPath+"output_ana_test.root","RECREATE");

  // data reader Object
  PndAnalysis* theAnalysis = new PndAnalysis();
  if (nevts==0) nevts = theAnalysis->GetEntries();
  
  //RhoCandLists for analysis
  RhoCandList piplus, piminus, lambda0, antiLambda0, proton, antiProton, xiplus, ximinus, xiSys;
  RhoCandList NotCombinedPiMinus, CombinedPiMinus, CombinedPiPlus, NotCombinedPiPlus;
  RhoCandList SelectedProton, SelectedAntiProton, SelectedPiMinus, SelectedPiPlus;
  RhoCandList Lambda0Fit, AntiLambda0Fit, XiMinusFit, XiPlusFit;
  RhoCandList mclist, all;

  //Dummy RhoCandidate
  RhoCandidate * dummyCand = new RhoCandidate();


  //***Mass selector
  double m0_lambda0= TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();
  cout<<"Mass of Lambda0: "<<m0_lambda0<<endl;
  RhoMassParticleSelector * lambdaMassSelector = new RhoMassParticleSelector("lambda0", m0_lambda0, 0.3);
 
  double m0_Xi = TDatabasePDG::Instance()->GetParticle("Xi-")->Mass();
  cout<<"Mass of Xi-: "<<m0_Xi<<endl;
  RhoMassParticleSelector * xiMassSelector = new RhoMassParticleSelector("Xi-", m0_Xi, 0.3);

  double m0_pbarpsystem = TDatabasePDG::Instance()->GetParticle("pbarpSystem")->Mass();

  
  double pbarmom = 2.7;
  double p_m0 = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
  TLorentzVector ini (0,0, pbarmom, sqrt(p_m0*p_m0+ pbarmom*pbarmom)+p_m0);
  TVector3 beamBoost = ini.BoostVector();
  
  PndRhoTupleQA qa(theAnalysis, pbarmom);
 
  int evt=-1;
  int index=0;

  while (theAnalysis->GetEvent() && ++evt<nevts){

    if ((evt%100)==0) cout << "evt "<< evt <<endl;
    cout << "Running event " << evt << endl;
		
    //***get MC list and store info
    theAnalysis->FillList(mclist, "McTruth");
    qa.qaMcList("", mclist, ntpMC);
    ntpMC->DumpData();
   	
		
//.........这里部分代码省略.........

{
  gSystem.Load("libPhysics");
  TLorentzVector lab(0.17,0.32,4.84,4.938);
  TLorentzVector zero(0.0,0.0,0.0,0.0);

  cout << "LAB = "<< lab[0]<< ", "<< lab[1]<< ", " << lab[2]<< ", " << lab[3]
       << ", rho = " << lab.Rho()<< ", theta = "<< lab.Theta()
       << ", phi = "<< lab.Phi()<< endl;

  TLorentzVector W = lab+zero;
  lab.Boost(-W.BoostVector());

  cout << "CM = "<< lab[0]<< ", "<< lab[1]<< ", " << lab[2]<< ", " << lab[3]
       << ", rho = " << lab.Rho()<< ", theta = "<< lab.Theta()
       << ", phi = "<< lab.Phi()
       << endl;

  /*

  Double_t masses[3]={0.938,0.135,0.135};

  TGenPhaseSpace event;
  event.SetDecay(W,3,masses);

  event.Generate();
  TLorentzVector *pProton = event.GetDecay(0);
  TLorentzVector *pPi0 = event.GetDecay(1);
  TLorentzVector *pPi02 = event.GetDecay(2);
  
  TH1D *h = new TH1D("his","theta",100,0,180);
  h->Fill(pProton->Theta()*57.3);

void angular::Loop()
{
   if (fChain == 0) return;

   TH1F *Hin_cosThetaStarLead = new TH1F("Hin_cosThetaStarLead","Hin_cosThetaStarLead",100,-1,1);
   TH1F *Hin_cosThetaStarSublead = new TH1F("Hin_cosThetaStarSublead","Hin_cosThetaStarSublead",100,-1,1);
   TH1F *Hin_deltaEta = new TH1F("Hin_deltaEta","Hin_deltaEta",100,-5,5);
   TH1F *Hin_deltaPhi = new TH1F("Hin_deltaPhi","Hin_deltaPhi",50,-1,1);
   
   TH1F *Hout_cosThetaStarLead = new TH1F("Hout_cosThetaStarLead","Hout_cosThetaStarLead",100,-1,1);
   TH1F *Hout_cosThetaStarSublead = new TH1F("Hout_cosThetaStarSublead","Hout_cosThetaStarSublead",100,-1,1);
   TH1F *Hout_deltaEta = new TH1F("Hout_deltaEta","Hout_deltaEta",100,-5,5);
   TH1F *Hout_deltaPhi = new TH1F("Hout_deltaPhi","Hout_deltaPhi",50,-1,1);

   Long64_t nentries = fChain->GetEntriesFast();

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

      TLorentzVector myLeadGamma(0,0,0,0);  
      myLeadGamma.SetPtEtaPhiM(leadPt, leadEta, leadPhi, 0.);
      TLorentzVector mySubleadGamma(0,0,0,0);  
      mySubleadGamma.SetPtEtaPhiM(subleadPt, subleadEta, subleadPhi, 0.);

      TLorentzVector grav = myLeadGamma + mySubleadGamma;

      TLorentzVector lead_Gstar(myLeadGamma);
      lead_Gstar.Boost(-grav.BoostVector());
      TLorentzVector sublead_Gstar(mySubleadGamma);
      sublead_Gstar.Boost(-grav.BoostVector());

      float leadCosThStar = lead_Gstar.CosTheta();
      float subleadCosThStar = sublead_Gstar.CosTheta();

      float deltaEta = leadEta-subleadEta;
      float deltaPhi = myLeadGamma.DeltaPhi(mySubleadGamma);
      float cosDPhi = cos(deltaPhi);

      // for data
      if (puweight==0) puweight=1;
      
      if (fabs(mass-755)<20){ 
	Hin_cosThetaStarLead    -> Fill(leadCosThStar,puweight);
	Hin_cosThetaStarSublead -> Fill(subleadCosThStar,puweight);
	Hin_deltaEta->Fill(deltaEta,puweight);
	Hin_deltaPhi->Fill(cosDPhi,puweight);
      } else {
	Hout_cosThetaStarLead    -> Fill(leadCosThStar,puweight);
	Hout_cosThetaStarSublead -> Fill(subleadCosThStar,puweight);
	Hout_deltaEta->Fill(deltaEta,puweight);
	Hout_deltaPhi->Fill(cosDPhi,puweight);
      }
   }

   TFile outfile("outfile.root","RECREATE");
   Hout_cosThetaStarLead->Write();
   Hout_cosThetaStarSublead->Write();
   Hout_deltaEta->Write();
   Hout_deltaPhi->Write();
   Hin_cosThetaStarLead->Write();
   Hin_cosThetaStarSublead->Write();
   Hin_deltaEta->Write();
   Hin_deltaPhi->Write();
}

//.........这里部分代码省略.........
			h_lep_min2_E	-> Fill(lep_min2.E());
			h_lep_min2_Pt	-> Fill(ptlepm2);
			h_lep_min2_eta	-> Fill(lep_min2.Eta());
			h_lep_min2_phi	-> Fill(lep_min2.Phi());	

		//reconstructing the two lepton pairs Lorentz vectors
		lpair1 = lep_plus1 + lep_min1;
		lpair2 = lep_plus2 + lep_min2;

		//constructing 3-vectors in the lab frame
		lep_plus1_lab 	= lep_plus1.Vect();
		lep_plus2_lab 	= lep_plus2.Vect();	//.Vect() gives 3 vector from 4vector
		lep_min1_lab 	= lep_min1.Vect();	
		lep_min2_lab 	= lep_min2.Vect();
		lpair1_lab 	= lep_plus1_lab.Cross(lep_min1_lab);	
		lpair2_lab 	= lep_plus2_lab.Cross(lep_min2_lab);

// 		cout << " pt of lepton pair1 on rest frame is: "<< lpair1.Perp()<<endl;


	   	//Filling up Histograms with angles defined in the lab frame
		h_angle_lab_pair1 -> Fill(lep_plus1_lab.Angle(lep_min1_lab));
		h_angle_lab_pair2 -> Fill(lep_plus2_lab.Angle(lep_min2_lab));


       		//Filling up histograms with variables from articles
       		h_angle_lab_deleta1	-> Fill(fabs(lep_min1.Eta() - lep_plus1.Eta()));
       		h_angle_lab_delphi1	-> Fill(fabs(lep_min1.DeltaPhi(lep_plus1)));
       		h_angle_lab_deleta2	-> Fill(fabs(lep_min2.Eta() - lep_plus2.Eta()));
       		h_angle_lab_delphi2	-> Fill(fabs(lep_min2.DeltaPhi(lep_plus2)));


	   	//Looking at the Higgs rest frame
	   	TVector3 boost_rH 	= -rec_H.BoostVector(); //NOTE the minus sign! WHY - sign???
	   	TVector3 boost_rZ1	= -Z1.BoostVector();
	   	TVector3 boost_rZ2	= -Z2.BoostVector();	
	   	Higgs_rest	= rec_H;
	   	Z1_rH		= Z1;
	   	Z2_rH		= Z2;
	   	lep_p1_rH	= lep_plus1;	//
	   	lep_m1_rH	= lep_min1;
	   	lep_p2_rH	= lep_plus2;
	   	lep_m2_rH	= lep_min2;
	   	lep_p1_rZ1	= lep_plus1;
	   	lep_m2_rZ2	= lep_min2;
	   	lep_p2_rZ2	= lep_plus2;
	   	lep_m1_rZ1	= lep_min1;

	   	//Boosting vectors to the Higgs rest frame
	   	Higgs_rest.Boost(boost_rH);
	   	Z1_rH.Boost(boost_rH);
	   	Z2_rH.Boost(boost_rH);
	   	lep_p1_rH.Boost(boost_rH);
	   	lep_m1_rH.Boost(boost_rH);
	   	lep_p2_rH.Boost(boost_rH);
	   	lep_m2_rH.Boost(boost_rH);

	   	//Boosting leptons to Z rest frames
	   	lep_p1_rZ1.Boost(boost_rZ1);
	   	lep_m1_rZ1.Boost(boost_rZ1);
	   	lep_p2_rZ2.Boost(boost_rZ2);
	   	lep_m2_rZ2.Boost(boost_rZ2);

	   	//Setting 3Vectors in Higgs rest frame
	   	Z3_1_rH		= Z1_rH.Vect();
	   	Z3_2_rH		= Z2_rH.Vect();

//.........这里部分代码省略.........
    tree->SetBranchAddress("pGenAntiNu",      &pGenAntiNu);
 
    int nEvents = (int)tree->GetEntries();
    //nEvents = 5000;
    int EventCounter = 0;
    cout << "Anzahl Ereignisse: " << nEvents << endl;
    for(int iEvent=1; iEvent<nEvents;iEvent++){

        tree->GetEntry(iEvent);

        EventCounter++;

  
        if(iEvent%10000 == 1)
        {
            cout << "Event " << iEvent << endl;
        }
        
        EventIsGood = 0;
        
// GENERATOR THETA

        w_A = 0;
        w_N = 0;
        
        *pTTbar=(*pTop+*pAntiTop);
        *pTopBoosted = *pTop;
        *pAntiTopBoosted = *pAntiTop;
        *pLeptonPlusBoosted = *pLeptonPlus;
        *pLeptonMinusBoosted = *pLeptonMinus;
        *pBBoosted = *pBQuark;
        *pBbarBoosted = *pBbarQuark;

        pAntiTopBoosted->Boost(-pTTbar->BoostVector());
        pTopBoosted->Boost(-pTTbar->BoostVector());
        pLeptonPlusBoosted->Boost(-pTop->BoostVector());
        pLeptonMinusBoosted->Boost(-pAntiTop->BoostVector());
        CosThetaPlus = cos(pLeptonPlusBoosted->Angle(pTopBoosted->Vect()));
        CosThetaMinus = cos(pLeptonMinusBoosted->Angle(pAntiTopBoosted->Vect()));

        pBBoosted->Boost(-pTop->BoostVector());
        pBbarBoosted->Boost(-pAntiTop->BoostVector());

        CosLeptonAngleD = cos(pLeptonPlusBoosted->Angle(pLeptonMinusBoosted->Vect()));

        double Nenner = 1 - 0.256*CosThetaPlus*CosThetaMinus;
        
        w_A = (-CosThetaPlus*CosThetaMinus)/Nenner;
        w_N = 1./Nenner;
        
        w_LL = (1-CosThetaPlus*CosThetaMinus-CosThetaPlus+CosThetaMinus)/Nenner;
        w_LR = (1+CosThetaPlus*CosThetaMinus-CosThetaPlus-CosThetaMinus)/Nenner;
        w_RR = (1-CosThetaPlus*CosThetaMinus+CosThetaPlus-CosThetaMinus)/Nenner;
        w_RL = (1+CosThetaPlus*CosThetaMinus+CosThetaPlus+CosThetaMinus)/Nenner;
        
        histogram__gen_A->Fill(CosThetaPlus, CosThetaMinus, w_A);
        histogram__gen_N->Fill(CosThetaPlus, CosThetaMinus, w_N);
        
        histogram__gen_LL->Fill(CosThetaPlus, CosThetaMinus, w_LL);
        histogram__gen_LR->Fill(CosThetaPlus, CosThetaMinus, w_LR);
        histogram__gen_RR->Fill(CosThetaPlus, CosThetaMinus, w_RR);
        histogram__gen_RL->Fill(CosThetaPlus, CosThetaMinus, w_RL);
        
        histogram__gen_Correlation->Fill(CosThetaPlus, CosThetaMinus);

        if(numberOfLeptons == 2) 

void number_of_particles_leaving_GEM_hits_boxgen(TString pre="", int nevts=0, double mom=4.1){
	  TDatabasePDG::Instance()-> AddParticle("pbarpSystem","pbarpSystem", 1.9, kFALSE, 0.1, 0,"", 88888);

	  TStopwatch timer;

	  if (pre==""){
		  //Output File
		  TString OutputFile = "test_analysis_output.root";
		  TString outPath = "";
		  //Input simulation Files
		 TString inPIDFile = "pid_complete.root";
		 TString inParFile = "simparams.root";
	  }
	  else {
		  //Output File
		  TString outPath = pre + "_";
		  TString OutputFile = pre + "_test_analysis_output.root";

		  //Input simulation Files
		  TString inPIDFile = pre + "_pid_complete.root";
		  TString inParFile = pre + "_simparams.root";
	  }

	  TString PIDParFile = TString( gSystem->Getenv("VMCWORKDIR")) + "/macro/params/all.par";

	  //Initialization
	  FairLogger::GetLogger()->SetLogToFile(kFALSE);
	  FairRunAna* RunAna = new FairRunAna();
	  FairRuntimeDb* rtdb = RunAna->GetRuntimeDb();
	  RunAna->SetInputFile(inPIDFile);

	  //setup parameter database
	  FairParRootFileIo* parIo = new FairParRootFileIo();
	  parIo->open(inParFile);
	  FairParAsciiFileIo* parIoPID = new FairParAsciiFileIo();
	  parIoPID->open(PIDParFile.Data(),"in");

	  rtdb->setFirstInput(parIo);
	  rtdb->setSecondInput(parIoPID);
	  rtdb->setOutput(parIo);

	  RunAna->SetOutputFile(OutputFile);
	  RunAna->Init();


	  /*************************************************************************
	   *  Create new ntuple and fill them with information
	   ************************************************************************/

	  //*** create tuples
	  RhoTuple * ntpPiMinus = new RhoTuple("ntpPiMinus", "PiMinus info");
	  RhoTuple * ntpPiPlus = new RhoTuple("ntpPiPlus", "PiPlus info");
	  RhoTuple * ntpKaonMinus = new RhoTuple("ntpKaonMinus", "KaonMinus info");
	  RhoTuple * ntpKaonPlus = new RhoTuple("ntpKaonPlus", "KaonPlus info");
	  RhoTuple * ntpProton = new RhoTuple("ntpProton", "Proton info");
	  RhoTuple * ntpAntiProton = new RhoTuple("ntpAntiProton", "Antiproton info");


	  //Create output file
	  TFile *out = TFile::Open(outPath+"test_output_ana.root","RECREATE");

	  // data reader Object
	  PndAnalysis* theAnalysis = new PndAnalysis();
	  if (nevts==0) nevts = theAnalysis->GetEntries();


	  //RhoCandLists for analysis
	  RhoCandList piplus, piminus, proton, antiproton, kaonminus, kaonplus;

	  RhoCandidate * dummyCand = new RhoCandidate(); //dummy candidate for empty candidate usage


	  double p_m0 = TDatabasePDG::Instance()->GetParticle("proton")->Mass();


	  TLorentzVector ini (0,0, mom, sqrt(p_m0*p_m0+ mom*mom)+p_m0);
	  TVector3 beamBoost = ini.BoostVector();

	  PndRhoTupleQA qa(theAnalysis, mom);

	  int evt=-1;
	  while (theAnalysis->GetEvent() && ++evt<nevts){

	    if ((evt%100)==0) cout << "evt "<< evt <<endl;


	    TString PidSelection = "PidAlgoIdealCharged";//"PidAlgoMvd;PidAlgoStt;PidAlgoDrc"; to change from ideal PID to realistic PID uncomment this!


	    //***Selection with no PID info
	    theAnalysis->FillList(piminus, "PionBestMinus", PidSelection);
	    theAnalysis->FillList(piplus, "PionBestPlus", PidSelection);
	    theAnalysis->FillList(kaonminus, "KaonBestMinus", PidSelection);
	    theAnalysis->FillList(kaonplus, "KaonBestPlus", PidSelection);
	    theAnalysis->FillList(proton, "ProtonBestPlus", PidSelection);
	    theAnalysis->FillList(antiproton, "ProtonBestMinus", PidSelection);


	    //Get piminus information
	    ntpPiMinus->Column("ev",     (Float_t) evt);
//.........这里部分代码省略.........

// Get corrections from table
float qq2vvEWKcorrections::getqq2WWEWKCorr(
                       float ptl1 , float etal1 , float phil1 , float idl1 , // lepton from 1st W
                       float ptl2 , float etal2 , float phil2 , float idl2 , // lepton from 2nd W
                       float ptv1 , float etav1 , float phiv1 ,             // neutrino from 1st W
                       float ptv2 , float etav2 , float phiv2 ,             // neutrino from 2nd W
                       float x1   , float x2 ,                              // parton x-Bjorken
                       float id1  , float id2                               // parton PDG id's
//                       int   TypeFirst = 1 ,                                // ??????????????????????????
//                        float Energy = 8000.
                     ) {

 float Energy = 8000.;
 // Create lepton and neutrino vectors
 TLorentzVector l1;
 TLorentzVector l2;
 TLorentzVector v1;
 TLorentzVector v2;

 float me   = 0.001*0.510998928 ;
 float mmu  = 0.001*105.6583715 ;
 float mtau = 0.001*1776.82     ;

 float ml1(0) , ml2(0) ;

 if      ( idl1 == 11 ) ml1 = me   ;
 else if ( idl1 == 13 ) ml1 = mmu  ;
 else if ( idl1 == 15 ) ml1 = mtau ;
 if      ( idl2 == 11 ) ml2 = me   ;
 else if ( idl2 == 13 ) ml2 = mmu  ;
 else if ( idl2 == 15 ) ml2 = mtau ;

 l1.SetPtEtaPhiM(ptl1,etal1,phil1,ml1);
 l2.SetPtEtaPhiM(ptl2,etal2,phil2,ml2);
 v1.SetPtEtaPhiM(ptv1,etav1,phiv1,0.0);
 v2.SetPtEtaPhiM(ptv2,etav2,phiv2,0.0);

 TLorentzVector W1 = l1+v1; // W1
 TLorentzVector W2 = l2+v2; // W2
 TLorentzVector WW = W1+W2;

 //---- FIX
 //---- swap 1 <->2 for leptons and neutrinos, to get correctly association l-v:
 //---- use as test invariatn mass of di-leptons
 //----    needed because leptons are ordered by "pt" and not by "mother" in the input! 
 if ((W1.M() - 80.385) > 0.1 && (W2.M() - 80.385) > 0.1) {
  W1 = l1+v2; // W1
  W2 = l2+v1; // W2
  WW = W1+W2;
 }
 //---- end FIX


 float M_12   = 80.385 , M_22 = 80.385 ;

 TLorentzVector p1; p1.SetPxPyPzE(0.,0.,Energy*x1,Energy*x1);
 TLorentzVector p2; p2.SetPxPyPzE(0.,0.,-Energy*x2,Energy*x2);

 //S-HAT
 float s_hat = pow(WW.E(),2)-pow(WW.Px(),2)-pow(WW.Py(),2)-pow(WW.Pz(),2); // ScalarProd(WW) (p1+p2)^2 = (p3+p4)^2 ~ +2*p1*p2
  //T_HAT
  //float t_hat2 = TypeFirst == 1 ? p1.M()*p1.M() + M_12*M_12 - 2*( p1.E()*W1.E() - p1.Px()*W1.Px() - p1.Py()*W1.Py() - p1.Pz()*W1.Pz() ) :
  //                                p1.M()*p1.M() + M_22*M_22 - 2*( p1.E()*W2.E() - p1.Px()*W2.Px() - p1.Py()*W2.Py() - p1.Pz()*W2.Pz() )   ; //T_HAT LO
 float la1 = sqrt( pow(s_hat,2) );            //la = sqrt( pow(a,2)+pow(b,2)+pow(c,2)-2*(a*b+a*c+b*c) );
 float la2 = sqrt( pow(s_hat,2) + pow(M_12,2) + pow(M_22,2) - 2*(s_hat*M_12 + s_hat*M_22 + M_12*M_22) );
  //  Boost: boost ext. momenta in CM frame of W1,W2
 TLorentzVector W1_b = W1, W2_b = W2;
 TLorentzVector p1_b = p1, p2_b = p2;
 W1_b.Boost( -WW.BoostVector() );
 W2_b.Boost( -WW.BoostVector() );
 p1_b.Boost( -WW.BoostVector() );
 p2_b.Boost( -WW.BoostVector() );

 //  Uni-vector
 TLorentzVector ee1 = p1_b*(1./sqrt( pow(p1_b.X(),2)+pow(p1_b.Y(),2)+pow(p1_b.Z(),2) ));
 TLorentzVector ee2 = p2_b*(1./sqrt( pow(p2_b.X(),2)+pow(p2_b.Y(),2)+pow(p2_b.Z(),2) ));
 TLorentzVector z1  = W1_b*(1./sqrt( pow(W1_b.X(),2)+pow(W1_b.Y(),2)+pow(W1_b.Z(),2) ));
 TLorentzVector z2  = W2_b*(1./sqrt( pow(W2_b.X(),2)+pow(W2_b.Y(),2)+pow(W2_b.Z(),2) ));
 //  "effective" beam axis
 float abse = sqrt( pow(ee1.X()-ee2.X(),2) + pow(ee1.Y()-ee2.Y(),2) + pow(ee1.Z()-ee2.Z(),2) );
 TLorentzVector ee = (ee1-ee2) * (1. / abse);
  //  "effective" scattering angle
 float costh = ee.X()*z1.X()+ee.Y()*z1.Y()+ee.Z()*z1.Z();
  //  final T_HAT
 float t_hat= M_12 - (1./2.)*(s_hat+M_12-M_22) + (1/(2*s_hat))*la1*la2*costh; //Mz-1/2*s+1/(2*s)*la(s,0,0)*la(s,mZ,mZ)*costh  or: (p1-p3)^2 = (p2-p4)^2 ~ -2*p1*p3
  //Quark Type
 int quarkType = -1.;
 if(      fabs(id1)==2 && fabs(id2)==2 ) quarkType=0; //delta_uub
 else if( fabs(id1)==1 && fabs(id2)==1 ) quarkType=1; //delta_ddb
 else if( fabs(id1)==5 && fabs(id2)==5 ) quarkType=2; //delta_bbb
 else if( fabs(id1)==4 && fabs(id2)==4 ) quarkType=0; // cc as delta_buu
 else if( fabs(id1)==3 && fabs(id2)==3 ) quarkType=1; // ss as delta_bdd

 // Extracting corrections
 float EWK_w = 1. ;
//  std::cout << " quarkType = " << quarkType << " id1 = " << id1 << " id2 = " << id2 << std::endl;
 if( quarkType!=-1 ){
  float sqrt_s_hat = sqrt(s_hat);
  std::vector<std::vector<float> > EWK_w2_vec = findCorrection( sqrt_s_hat, t_hat );
  float EWK_w2 = 1. + EWK_w2_vec[0][quarkType];
//.........这里部分代码省略.........

inline TVector3 kinematics::systemBoostVector( TLorentzVector* pa, TLorentzVector* pb )
{
	TLorentzVector pTmp = (*pa)+(*pb);
	return pTmp.BoostVector();
}

void calculateAngles(TLorentzVector thep4H, TLorentzVector thep4Z1, TLorentzVector thep4Lep11, TLorentzVector thep4Lep12, TLorentzVector thep4Z2, TLorentzVector thep4Lep21, TLorentzVector thep4Lep22, double& costheta1, double& costheta2, double& phi, double& costhetastar, double& phistar1){
	
	
	//std::cout << "In calculate angles..." << std::endl;
	
	double norm;
	
	TVector3 boostX = -(thep4H.BoostVector());
	TLorentzVector thep4Z1inXFrame( thep4Z1 );
	TLorentzVector thep4Z2inXFrame( thep4Z2 );	
	thep4Z1inXFrame.Boost( boostX );
	thep4Z2inXFrame.Boost( boostX );
	TVector3 theZ1X_p3 = TVector3( thep4Z1inXFrame.X(), thep4Z1inXFrame.Y(), thep4Z1inXFrame.Z() );
	TVector3 theZ2X_p3 = TVector3( thep4Z2inXFrame.X(), thep4Z2inXFrame.Y(), thep4Z2inXFrame.Z() );
	
	///////////////////////////////////////////////
	// check for z1/z2 convention, redefine all 4 vectors with convention
	///////////////////////////////////////////////	
	TLorentzVector p4H, p4Z1, p4M11, p4M12, p4Z2, p4M21, p4M22;

	p4Z1 = thep4Z1; p4M11 = thep4Lep11; p4M12 = thep4Lep12;
	p4Z2 = thep4Z2; p4M21 = thep4Lep21; p4M22 = thep4Lep22;
	costhetastar = theZ1X_p3.CosTheta();
	
	// now helicity angles................................
	// ...................................................
	TVector3 boostZ1 = -(p4Z1.BoostVector());
	TLorentzVector p4Z2Z1(p4Z2);
	p4Z2Z1.Boost(boostZ1);
	//find the decay axis
	/////TVector3 unitx_1 = -Hep3Vector(p4Z2Z1);
	TVector3 unitx_1( -p4Z2Z1.X(), -p4Z2Z1.Y(), -p4Z2Z1.Z() );
	norm = 1/(unitx_1.Mag());
	unitx_1*=norm;
	//boost daughters of z2
	TLorentzVector p4M21Z1(p4M21);
	TLorentzVector p4M22Z1(p4M22);
	p4M21Z1.Boost(boostZ1);
	p4M22Z1.Boost(boostZ1);
	//create z and y axes
	/////TVector3 unitz_1 = Hep3Vector(p4M21Z1).cross(Hep3Vector(p4M22Z1));
	TVector3 p4M21Z1_p3( p4M21Z1.X(), p4M21Z1.Y(), p4M21Z1.Z() );
	TVector3 p4M22Z1_p3( p4M22Z1.X(), p4M22Z1.Y(), p4M22Z1.Z() );
	TVector3 unitz_1 = p4M21Z1_p3.Cross( p4M22Z1_p3 );
	norm = 1/(unitz_1.Mag());
	unitz_1 *= norm;
	TVector3 unity_1 = unitz_1.Cross(unitx_1);
	
	//caculate theta1
	TLorentzVector p4M11Z1(p4M11);
	p4M11Z1.Boost(boostZ1);
	TVector3 p3M11( p4M11Z1.X(), p4M11Z1.Y(), p4M11Z1.Z() );
	TVector3 unitM11 = p3M11.Unit();
	double x_m11 = unitM11.Dot(unitx_1); double y_m11 = unitM11.Dot(unity_1); double z_m11 = unitM11.Dot(unitz_1);
	TVector3 M11_Z1frame(y_m11, z_m11, x_m11);
	costheta1 = M11_Z1frame.CosTheta();
	//std::cout << "theta1: " << M11_Z1frame.Theta() << std::endl;
	//////-----------------------old way of calculating phi---------------/////////
	phi = M11_Z1frame.Phi();
	
	//set axes for other system
	TVector3 boostZ2 = -(p4Z2.BoostVector());
	TLorentzVector p4Z1Z2(p4Z1);
	p4Z1Z2.Boost(boostZ2);
	TVector3 unitx_2( -p4Z1Z2.X(), -p4Z1Z2.Y(), -p4Z1Z2.Z() );
	norm = 1/(unitx_2.Mag());
	unitx_2*=norm;
	//boost daughters of z2
	TLorentzVector p4M11Z2(p4M11);
	TLorentzVector p4M12Z2(p4M12);
	p4M11Z2.Boost(boostZ2);
	p4M12Z2.Boost(boostZ2);
	TVector3 p4M11Z2_p3( p4M11Z2.X(), p4M11Z2.Y(), p4M11Z2.Z() );
	TVector3 p4M12Z2_p3( p4M12Z2.X(), p4M12Z2.Y(), p4M12Z2.Z() );
	TVector3 unitz_2 = p4M11Z2_p3.Cross( p4M12Z2_p3 );
	norm = 1/(unitz_2.Mag());
	unitz_2*=norm;
	TVector3 unity_2 = unitz_2.Cross(unitx_2);
	//calcuate theta2
	TLorentzVector p4M21Z2(p4M21);
	p4M21Z2.Boost(boostZ2);
	TVector3 p3M21( p4M21Z2.X(), p4M21Z2.Y(), p4M21Z2.Z() );
	TVector3 unitM21 = p3M21.Unit();
	double x_m21 = unitM21.Dot(unitx_2); double y_m21 = unitM21.Dot(unity_2); double z_m21 = unitM21.Dot(unitz_2);
	TVector3 M21_Z2frame(y_m21, z_m21, x_m21);
	costheta2 = M21_Z2frame.CosTheta();
	
	// calculate phi
	//calculating phi_n
	TLorentzVector n_p4Z1inXFrame( p4Z1 );
	TLorentzVector n_p4M11inXFrame( p4M11 );
	n_p4Z1inXFrame.Boost( boostX );
	n_p4M11inXFrame.Boost( boostX );        
	TVector3 n_p4Z1inXFrame_unit = n_p4Z1inXFrame.Vect().Unit();
	TVector3 n_p4M11inXFrame_unit = n_p4M11inXFrame.Vect().Unit();  
	TVector3 n_unitz_1( n_p4Z1inXFrame_unit );
	//// y-axis is defined by neg lepton cross z-axis
	//// the subtle part is here...
	//////////TVector3 n_unity_1 = n_p4M11inXFrame_unit.Cross( n_unitz_1 );
	TVector3 n_unity_1 = n_unitz_1.Cross( n_p4M11inXFrame_unit );
//.........这里部分代码省略.........