C++ TVector2::SetMagPhi方法代码示例

void ttreco2( TString process ) {

	cout << "\n\n" << process << "\n" <<  endl; 

	TFile* fshape  = new TFile( "/afs/cern.ch/user/j/jgarciaf/mimick/mlb.root" ); 

	TH1F* shapemlb = (TH1F*) fshape->Get( "mlb" );

	MassReconstructor theMass( 100, shapemlb );  


	//TFile myfile( "../minitrees/" + inputdir + "/TTDM/" + process + ".root", "update" ); 
	TFile myfile( "../../../../public/minitrees_week-1_with-ttReco/" + process + ".root", "update" ); 
	//TFile myfile( "../minitrees/week-1-checks/" + process + ".root", "update" ); 

   	TTreeReader myreader( "latino", &myfile );

	TTree* mytree = (TTree*) myfile.Get( "latino" );


	//----- read -------------------------------------------------------

   	TTreeReaderValue<float> metPfType1   ( myreader, "metPfType1"    );
   	TTreeReaderValue<float> metPfType1Phi( myreader, "metPfType1Phi" );

   	TTreeReaderValue<float> lep1pt  ( myreader, "lep1pt"   );
   	TTreeReaderValue<float> lep1eta ( myreader, "lep1eta"  );
   	TTreeReaderValue<float> lep1phi ( myreader, "lep1phi"  );
   	TTreeReaderValue<float> lep1mass( myreader, "lep1mass" );

	TTreeReaderValue<float> lep2pt  ( myreader, "lep2pt"   );
   	TTreeReaderValue<float> lep2eta ( myreader, "lep2eta"  );
   	TTreeReaderValue<float> lep2phi ( myreader, "lep2phi"  );
   	TTreeReaderValue<float> lep2mass( myreader, "lep2mass" );

   	TTreeReaderValue<std::vector<float>> jet_pt ( myreader, "jet_pt"  );
   	TTreeReaderValue<std::vector<float>> jet_eta( myreader, "jet_eta" );
   	TTreeReaderValue<std::vector<float>> jet_phi( myreader, "jet_phi" );

   	TTreeReaderValue<std::vector<float>> bjet30csvv2m_pt ( myreader, "bjet30csvv2m_pt"  );
   	TTreeReaderValue<std::vector<float>> bjet30csvv2m_eta( myreader, "bjet30csvv2m_eta" );
   	TTreeReaderValue<std::vector<float>> bjet30csvv2m_phi( myreader, "bjet30csvv2m_phi" );

	//TTreeReaderValue<float> top1pt  ( myreader, "top1pt_gen"   );
   	//TTreeReaderValue<float> top1eta ( myreader, "top1eta_gen"  );
   	//TTreeReaderValue<float> top1phi ( myreader, "top1phi_gen"  );

	//TTreeReaderValue<float> top2pt  ( myreader, "top2pt_gen"   );
   	//TTreeReaderValue<float> top2eta ( myreader, "top2eta_gen"  );
   	//TTreeReaderValue<float> top2phi ( myreader, "top2phi_gen"  );


	//----- write ------------------------------------------------------

	float topRecoW;
	float darkpt  ; 

	TBranch* b_topRecoW = mytree -> Branch( "topRecoW", &topRecoW, "topRecoW/F" );
	TBranch* b_darkpt   = mytree -> Branch( "newdarkpt"  , &darkpt  , "darkpt/F"   );

	TBranch* b_mimick = mytree -> Branch( "mimick"  , &mimick  , "mimick/F"   );

	

	//----- loop -------------------------------------------------------


	int nentries = myreader.GetEntries(1); 

	//if ( nentries > 10 ) nentries = 10; 

	for ( Long64_t ievt = 0; ievt < nentries; ievt++ ) {

		if( ievt%10000 == 0 ) cout << "\n\n ievt: " << ievt << endl;

		myreader.SetEntry(ievt);

		//--- MET

		TVector2 MET;

		MET.SetMagPhi( *metPfType1, *metPfType1Phi );

		myreader.SetEntry(ievt);

		//--- leptons

		TLorentzVector l1, l2;   // top1gen, top2gen; 

		l1.SetPtEtaPhiM( *lep1pt, *lep1eta, *lep1phi, *lep1mass ); 
  		l2.SetPtEtaPhiM( *lep2pt, *lep2eta, *lep2phi, *lep2mass );

		//top1gen.SetPtEtaPhiM( *top1pt, *top1eta, *top1phi, 173.34 );
		//top2gen.SetPtEtaPhiM( *top2pt, *top2eta, *top2phi, 173.34 );


		//--- jets

		std::vector<TLorentzVector> jets;
 		std::vector<TLorentzVector> bjets;
//.........这里部分代码省略.........

void PMCSZCand::SetA(PMCSEMObj &elec1, PMCSEMObj &elec2) {

  // First we need to calculate the thrust axis
  //
  TVector2 e1(elec1.ppx(),elec1.ppy());
  TVector2 e2(elec2.ppx(),elec2.ppy());

  // Calculate the two phi angles
  double phi1=e1.Phi();
  double phi2=e2.Phi();

  // Order in phi
  if (phi1<phi2) {
    TVector2 dummy=e1;
    e1=e2;
    e2=dummy;
    double dummy2=phi1;
    phi1=phi2;
    phi2=dummy2;
  }

  // Calculate lengths as well
  double len1=e1.Mod();
  double len2=e2.Mod();

  // Good old Newton
  //
  // initial guess
  TVector2 bisector=(e1.Unit()+e2.Unit()).Unit();
  double alpha=bisector.Phi()-TMath::Pi()/2.;
  double alphaBackup=alpha;
  int nIt=0;
  // iterate
  double oldAlpha=9999.;
  while (fabs(alpha-oldAlpha)>0.000001) {
    oldAlpha=alpha;
    double f=len2*sin(phi2-oldAlpha)+len1*sin(oldAlpha-phi1);
    double fp=-len2*cos(phi2-oldAlpha)+len1*cos(oldAlpha-phi1);
    alpha=oldAlpha-f/fp;
    nIt++;
    if (nIt>1000) {
      cout<<"Newton did not converge in search for thrust axis"<<endl;
      alpha=alphaBackup;
      break;
    }
  }

  // Build unit vector
  TVector2 r;
  r.SetMagPhi(1.,alpha);
  TVector2 rPerp=r.Rotate(TMath::Pi()/2.);

  // Checks
  if (fabs(e1*rPerp-e2*rPerp)>0.001) {
    cout<<"Looks like Newton was imprecise in calculation of thrust axis: "<<e1*rPerp<<" "<<e2*rPerp<<endl;
  }

  // now we can do the projections
  TVector2 Z=e1+e2;
  _pat=Z*rPerp;
  _pal=Z*r;
}

Bool_t test::Process(Long64_t entry)
{
  // The Process() function is called for each entry in the tree (or possibly
  // keyed object in the case of PROOF) to be processed. The entry argument
  // specifies which entry in the currently loaded tree is to be processed.
  // It can be passed to either test::GetEntry() or TBranch::GetEntry()
  // to read either all or the required parts of the data. When processing
  // keyed objects with PROOF, the object is already loaded and is available
  // via the fObject pointer.
  //
  // This function should contain the "body" of the analysis. It can contain
  // simple or elaborate selection criteria, run algorithms on the data
  // of the event and typically fill histograms.
  //
  // The processing can be stopped by calling Abort().
  //
  // Use fStatus to set the return value of TTree::Process().
  //
  // The return value is currently not used.
  if (entry % 1000 == 0){
    std::cout<<"Entry: "<<entry<<std::endl;
  }

  GetEntry(entry);

  w=WT00;

  TLorentzVector * nVec = new TLorentzVector() ;
  TLorentzVector * eVec = new TLorentzVector() ;
  TLorentzVector * wVec = new TLorentzVector() ;


  nVec->SetPxPyPzE(Px_d2, Py_d2, Pz_d2, E_d2);
  eVec->SetPxPyPzE(Px_d1, Py_d1, Pz_d1, E_d1);
  wVec->SetPxPyPzE(Px_V, Py_V, Pz_V, E_V);

  if (eVec->Pt()>lep_ET){
    if (fabs(eVec->Eta())<lep_eta){
      if (fabs(eVec->Eta())<lep_crack_1 || fabs(eVec->Eta())>lep_crack_2){
        //std::cout<<"met: "<<met<<std::endl;
        if (bZ){
          if (fabs(nVec->Eta())<lep_crack_1 || fabs(nVec->Eta())>lep_crack_2){
            TVector2 * MET = new TVector2(0.,0.);

            TVector2 sumLepton(Px_d1+Px_d2, Py_d1+Py_d2);

            SetTreeVariables(wVec, wVec, &sumLepton, MET, MET, w, 0, iTree);

            delete MET;
          }
        }
        else{
          double met = nVec->Pt();
          double met_phi = nVec->Phi();
          TVector2 * MET = new TVector2(0.,0.);

          TVector2 sumLepton(Px_d1, Py_d1);

          MET->SetMagPhi(met,met_phi);
          SetTreeVariables(wVec, wVec, &sumLepton, MET, MET, w, 0, iTree);

          delete MET;
        }
      }
    }
  }
  delete nVec;
  delete eVec;
  delete wVec;
  return kTRUE;

}

void fillCategory(EventTree* tree, TH1F* hist, double weight){
	int EleP = 0;
	int EleM = 0;
	int MuP = 0;
	int MuM = 0;
	int TauP = 0;
	int TauM = 0;
	int ElePfid = 0;
	int EleMfid = 0;
	int MuPfid = 0;
	int MuMfid = 0;
	int nNufid = 0;

	int nPhofid = 0;
	int nJetfid = 0;
	int nBJetfid = 0;

	TVector2 MET = TVector2(0,0);
	TVector2 tempNu = TVector2(0,0);

	for( int mcI = 0; mcI < tree->nMC_; ++mcI){
	  if(abs(tree->mcMomPID->at(mcI))==24 && tree->mcParentage->at(mcI)==10){
	    if( tree->mcPID->at(mcI) == 11 ) EleP = 1;
	    if( tree->mcPID->at(mcI) == -11 ) EleM = 1;
	    if( tree->mcPID->at(mcI) == 13 ) MuP = 1;
	    if( tree->mcPID->at(mcI) == -13 ) MuM = 1;
	    if( tree->mcPID->at(mcI) == 15) TauP = 1;
	    if( tree->mcPID->at(mcI) == -15) TauM = 1;
	  }
	  if((abs(tree->mcMomPID->at(mcI))==24 && tree->mcParentage->at(mcI)==10) || (abs(tree->mcMomPID->at(mcI))==15 && tree->mcParentage->at(mcI)==26)){		  
	    if( tree->mcPID->at(mcI) == 11 ) {
	      if (tree->mcPt->at(mcI) > 35 && (fabs(tree->mcEta->at(mcI)) < 2.5 && !(fabs(tree->mcEta->at(mcI)) > 1.4442 && fabs(tree->mcEta->at(mcI))<1.566))) ElePfid += 1;
	    }
	    if( tree->mcPID->at(mcI) == -11 ) {
	      if (tree->mcPt->at(mcI) > 35 && (fabs(tree->mcEta->at(mcI)) < 2.5 && !(fabs(tree->mcEta->at(mcI)) > 1.4442 && fabs(tree->mcEta->at(mcI))<1.566))) EleMfid += 1;
	    }
	    if( tree->mcPID->at(mcI) == 13 ) {
	      if (tree->mcPt->at(mcI) > 26 && fabs(tree->mcEta->at(mcI)) < 2.1) MuPfid += 1;
	    }
	    if( tree->mcPID->at(mcI) == -13 ) {
	      if (tree->mcPt->at(mcI) > 26 && fabs(tree->mcEta->at(mcI)) < 2.1) MuMfid += 1;
	    }
	  }
	  if( fabs(tree->mcPID->at(mcI)) == 12 || fabs(tree->mcPID->at(mcI)) == 14 || fabs(tree->mcPID->at(mcI)) == 16 ) {
	    if (tree->mcPt->at(mcI) > 20){
	      nNufid += 1;
	      //cout << nNufid << "\t" << tree->run_ << tree->lumis_ << tree->event_ << "\t" << mcI << "\t" <<  tree->mcPt->at(mcI) << "\t" << tree->mcEta->at(mcI) << "\t" << tree->mcPhi->at(mcI) << endl;
	    }
	    tempNu.SetMagPhi(tree->mcPt->at(mcI),tree->mcPhi->at(mcI));
	    MET += tempNu;
	    //	    cout << nNufid << "\t" << tree->run_ << tree->lumis_ << tree->event_ << "\t" << mcI << "\t" << tree->mcPID->at(mcI) << "\t" <<  tree->mcPt->at(mcI) << "\t" << tree->mcEta->at(mcI) << "\t" << tree->mcPhi->at(mcI) << "\t" << MET.Mod() << endl;
	  }

	  if(tree->mcPID->at(mcI) == 22 && 
	     (tree->mcParentage->at(mcI)==2 || tree->mcParentage->at(mcI)==10 || tree->mcParentage->at(mcI)==26) && 
	     tree->mcPt->at(mcI) > 25 && 
	     fabs(tree->mcEta->at(mcI)) < 1.4442){
	    nPhofid += 1;
	  }

	  if( abs(tree->mcPID->at(mcI)) < 6 && (abs(tree->mcMomPID->at(mcI))==24 || abs(tree->mcMomPID->at(mcI))==6 ) && tree->mcPt->at(mcI) > 30 && abs(tree->mcEta->at(mcI))<2.4 ) {
	    nJetfid += 1;
	    if (abs(tree->mcMomPID->at(mcI))==6 && abs(tree->mcPID->at(mcI))==5){
	      nBJetfid += 1;
	    }
	  }
	}
	
	hist->Fill(1.0, weight); // Total
	int nEle = EleP + EleM;
	int nMu = MuP + MuM;
	int nTau = TauP + TauM;
	if( nEle + nMu + nTau == 0) hist->Fill(2.0, weight); // All Had
	if( nEle + nMu + nTau == 1) hist->Fill(3.0, weight); // Single Lepton
	if( nEle + nMu + nTau == 2) hist->Fill(4.0, weight); // Di Lepton

	int nElefid = ElePfid + EleMfid;
	int nMufid = MuPfid + MuMfid;

	if ( nEle==1 && nMu==0 && nTau==0 ) hist->Fill(6.0, weight);
	if ( nEle==0 && nMu==1 && nTau==0 ) hist->Fill(7.0, weight);
	if ( nEle==0 && nMu==0 && nTau==1 ) hist->Fill(8.0, weight);
	if ( nElefid==1 && nMufid==0) hist->Fill(9.0, weight); //Ejets final state (generated)
	if ( nElefid==0 && nMufid==1) hist->Fill(10.0, weight); //Mujets final state (generated)

	// if ( nElefid==1 && nMufid==0) hist->Fill(6.0, weight); //Ejets final state (generated)
	// if ( nElefid==2 && nMufid==0) hist->Fill(7.0, weight); //Ejets final state (generated)
	// if ( nElefid==0 && nMufid==1) hist->Fill(8.0, weight); //Mujets final state (generated)
	// if ( nElefid==0 && nMufid==2) hist->Fill(9.0, weight); //Ejets final state (generated)
	// if ( nElefid==1 && nMufid==1) hist->Fill(10.0, weight); //Ejets final state (generated)

	int nJetsfid = 0;
	int nBJetsfid = 0;
	
	if ((nElefid + nMufid)==1){
	  // for ( int jetI = 0; jetI < tree->genJetPt_->size(); jetI++){
	  //   if (tree->genJetPt_->at(jetI) >= 30 && fabs(tree->genJetEta_->at(jetI)) < 2.4) nJetsfid += 1;	    
	  // }
	  for ( int jetI = 0; jetI < tree->nJet_; jetI++){
	    //	    if (tree->jetGenPt_->at(jetI) >= 30 && fabs(tree->jetGenEta_->at(jetI)) < 2.4) nJetsfid += 1;
//.........这里部分代码省略.........