C++ TFile::Save方法代码示例

void dR(){
  gROOT->Reset();
  gStyle->SetOptStat(1);

  ifstream dataphi;
  ifstream dataeta;
  
  float etaRPC;
  float etaCSC;

  float phiRPC;
  float phiCSC;

  dataphi.open("phi.txt");
  dataeta.open("eta.txt");

  std::string name = "dR";

  float pi = 3.14159265;
  float dphi = 0;
  float deta = 0;

  TH1F * histo = new TH1F("histo",name.c_str(),1000,0,0.5);
  TH2F * histo2D = new TH2F("histo2D",name.c_str(),400,-0.2,0.2,400,-0.2,0.2);

  while(!dataphi.eof()){
    dataphi >>phiRPC>>phiCSC;
    dataeta >>etaRPC>>etaCSC;
    dphi = phiRPC-phiCSC;
    deta = etaRPC-etaCSC;
    histo->Fill(sqrt(dphi*dphi+deta*deta));
    histo2D->Fill(dphi,deta);
  }

  TFile * theFile = new TFile("dR.root","RECREATE");
  theFile->cd();

  TCanvas * Ca0 = new TCanvas("Ca0","Canvas",1200,800);
 /* histo->SetLineWidth(3);
  histo->Draw();
  histo->SetMaximum(histo->GetMaximumStored());
  histo->GetXaxis()->SetTitle(name.c_str());*/

  gStyle->SetPalette(1);

  histo2D->GetXaxis()->SetTitle("dPhi");
  histo2D->GetYaxis()->SetTitle("dEta");
  histo2D->Draw();
  histo2D->SetDrawOption("color");
  histo2D->Write();
  Ca0->SaveAs((name+".png").c_str());
  theFile->Save();
  theFile->Close();
 
}

void scatterPloteta(){
  gROOT->Reset();
  gStyle->SetOptStat(1);

  ifstream data;
  
  float valueX;
  float valueY;

  string name = "eta";
  data.open((name+".txt").c_str());
  
  TH2F * histo = new TH2F("histo",name.c_str(),1000,-2.4,2.4,100,-2.4,2.4); //eta
  TH1F * histo1D = new TH1F("histo1D",name.c_str(),1000,-2.4,2.4); //eta

//  float pi = 3.14159265;
//  TH2F * histo = new TH2F("histo",name.c_str(),1000,-1*pi,pi,1000,-1*pi,pi); //phi
//  TH1F * histo1D = new TH1F("histo1D",name.c_str(),1000,-1*pi,pi); //phi

  while(!data.eof()){
    data >>valueX>>valueY;
    cout<<valueX<<" "<<valueY<<endl;
    histo->Fill(valueX,valueY);
    histo1D->Fill(valueX-valueY);
  }
  
  TCanvas * Ca0 = new TCanvas("Ca0","Canvas",1200,800);
  histo->SetLineWidth(3);

  TFile * theFile = new TFile("eta.root","RECREATE");
  theFile->cd();

  histo->Draw();
  histo1D->DrawNormalized("same");
  histo->GetXaxis()->SetTitle(("RPC"+name).c_str());
  histo->GetYaxis()->SetTitle(("CSC"+name).c_str());
  Ca0->SaveAs((name+".png").c_str());

  histo1D->Write();
  histo->Write();

  theFile->Save();
  theFile->Close();

  
}

void DalitzplotXi1820_XiPlus_lambda_k(){


	//*** Data input
	TString inputFile = "/home/ikp1/puetz/panda/mysimulations/analysis/XiPlus_Lambda_K/output_ana.root";
	TFile * data = new TFile(inputFile, "READ");

	TString outPath = "/home/ikp1/puetz/panda/mysimulations/analysis/XiPlus_Lambda_K/plots";
	TFile * out = new TFile(outPath+"/root-files/Dalitzplot_MC.root", "RECREATE");


	TTree * sim = (TTree*) data->Get("ntpMC");
	int nevents = sim->GetEntriesFast();


	TH2D * dalitz_Xilk = new TH2D("dalitz_Xilk", "Dalitz plot for MC; m^{2}(#Lambda^{0},K^{-})/GeV^{2}/c^{4}; m^{2}(#bar{#Xi}, K^{-})/GeV^{2}/c^{4}", 150,2.5,3.8,150,3.2,4.8);
	gStyle->SetOptStat(0);

	TLorentzVector lXi, lk, lla, lXi1820;
	TLorentzVector PXiK, PlaK, PXil;

	for (int n=0; n<nevents; n++){

		sim->GetEntry(n);

		double Eaxi = sim->GetLeaf("e")->GetValue(xi);
		double Ek = sim->GetLeaf("e")->GetValue(kaon);
		double Ela = sim->GetLeaf("e")->GetValue(lambda0);

		double Pxaxi = sim->GetLeaf("px")->GetValue(xi);
		double Pxk = sim->GetLeaf("px")->GetValue(kaon);
		double Pxla = sim->GetLeaf("px")->GetValue(lambda0);

		double Pyaxi = sim->GetLeaf("py")->GetValue(xi);
		double Pyk = sim->GetLeaf("py")->GetValue(kaon);
		double Pyla = sim->GetLeaf("py")->GetValue(lambda0);

		double Pzaxi = sim->GetLeaf("pz")->GetValue(xi);
		double Pzk = sim->GetLeaf("pz")->GetValue(kaon);
		double Pzla = sim->GetLeaf("pz")->GetValue(lambda0);

		lXi.SetPxPyPzE(Pxaxi, Pyaxi, Pzaxi, Eaxi);
		lk.SetPxPyPzE(Pxk, Pyk, Pzk, Ek);
		lla.SetPxPyPzE(Pxla, Pyla, Pzla, Ela);


		PXiK = lXi + lk;
		PlaK = lla + lk;
		PXil = lXi + lla;


		dalitz_Xilk->Fill(PlaK.M2(),PXiK.M2());

	}

	setPandaStyle();

	out->cd();

	dalitz_Xilk->Write();

	out->Save();


	TCanvas * c = new TCanvas("c", "Dalitz plot PHSP model", 0,0,1500,1000);
	dalitz_Xilk->Draw("COLZ");


	PandaSmartLabel("L");

	//****write histograms
	c->Print(outPath+"/png-files/Dalitzplots_MC.png");
	c->Print(outPath+"/pdf-files/Dalitzplots_MC.pdf");


}

//.........这里部分代码省略.........
      ntpXiPlus->DumpData();
	 }

    AntiLambda0Fit.Cleanup();
//    BestCandAntiLambda0.Cleanup();
	BestVtxFitXiPlus.clear();
	BestMassFitXiPlus.clear();



    //******* Xi+ Xi- System*****************************

    xiSys.Combine(XiPlusFit, XiMinusFit);
    xiSys.SetType(88888);

    for (int syscand=0; syscand<xiSys.GetLength(); ++syscand){

		ntpXiSys->Column("ev",     (Float_t) evt);
		ntpXiSys->Column("cand",    (Float_t) j);
		ntpXiSys->Column("ncand",   (Float_t) ximinus.GetLength());
		ntpXiSys->Column("McTruthMatch", (bool) theAnalysis->McTruthMatch(xiSys[syscand]));

		RhoCandidate * mother = xiSys[syscand]->TheMother();
		int moth = (mother==0x0) ? 88888 : mother->PdgCode();
		ntpXiSys->Column("Mother", (Float_t) moth);

		qa.qaP4("XiSys_", xiSys[syscand]->P4(), ntpXiSys);
		qa.qaComp("XiSys_", xiSys[syscand], ntpXiSys);
		qa.qaPoca("XiSys_", xiSys[syscand], ntpXiSys);


		RhoCandidate *  truth = xiSys[syscand]->GetMcTruth();
		TLorentzVector l;

		if (truth != 0x0){
//			qa.qaComp("McTruth_", truth, ntpXiSys);
			qa.qaVtx("McTruth_", truth, ntpXiSys);
			l = truth->P4();
		}
		else{
//			qa.qaComp("McTruth_", dummyCand, ntpXiSys);
			qa.qaVtx("McTruth_", dummyCand, ntpXiSys);
		}
		qa.qaP4("McTruth_", l, ntpXiSys);


		//4C-Fitter

		PndKinFitter fitter4c (xiSys[syscand]);
		fitter4c.Add4MomConstraint(ini);
		fitter4c.Fit();

		RhoCandidate * xiSysFit4c = xiSys[syscand]->GetFit();

		qa.qaFitter("4CFit_", &fitter4c, ntpXiSys);
		qa.qaComp("4cFit_", xiSysFit4c, ntpXiSys);
		qa.qaVtx("4CFit_", xiSysFit4c, ntpXiSys);


		ntpXiSys->DumpData();


    }
    XiMinusFit.Cleanup();
    XiPlusFit.Cleanup();
  }



   

  //Write output
  out->cd();

  ntpMC -> GetInternalTree()->Write();
  ntpPiMinus ->GetInternalTree()->Write();
  ntpPiPlus->GetInternalTree()->Write();
  ntpProton->GetInternalTree()->Write();
  ntpAntiProton->GetInternalTree()->Write();
  ntpLambda0->GetInternalTree()->Write();
  ntpAntiLambda0->GetInternalTree()->Write();
  ntpXiMinus->GetInternalTree()->Write();
  ntpXiPlus->GetInternalTree()->Write();
  ntpXiSys->GetInternalTree()->Write();

  out->Save();
  

  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
  cout<<endl<<endl;
  cout<<"Macro finisched successfully."<<endl;
  cout<<"Realtime: "<<rtime<<" s, CPU time: "<<ctime<<" s"<<endl;
  cout<<endl;

 
  exit(0);
 
}

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


	    }

	    ntpKaonPlus->DumpData();

//		Get Proton information
	    ntpProton->Column("ev", (int) evt);

	    for (int j=0; j<proton.GetLength(); ++j){


			//information about the mother and MCTruth Candidate

			TLorentzVector l;
			float costheta = -999.;

			RhoCandidate * truth = proton[j]->GetMcTruth();
			RhoCandidate * mother;
			if (truth)  mother = truth->TheMother();

			int moth = (mother==0x0) ? 88888 : mother->PdgCode();
			ntpProton->Column("Mother", (Int_t) moth);

			bool truthmatch = theAnalysis->McTruthMatch(proton[j]);
			ntpProton->Column("MCTruthMatch", (bool) truthmatch);

			int gemhit = GemHits(proton[j]);

			int count = 0;
			if (moth==88888 && gemhit==1 && truthmatch==1) count=1;

			ntpProton->Column("GemHit", (int) count, 0);


	    }

	    ntpProton->DumpData();

//			Get Antiproton
	    ntpAntiProton->Column("ev", (int) evt);

	    for (int j=0; j<antiproton.GetLength(); ++j){


			//information about the mother and MCTruth Candidate

			TLorentzVector l;
			float costheta = -999.;

			RhoCandidate * truth = antiproton[j]->GetMcTruth();
			RhoCandidate * mother;
			if (truth)  mother = truth->TheMother();

			int moth = (mother==0x0) ? 88888 : mother->PdgCode();
			ntpAntiProton->Column("Mother", (Int_t) moth);

			bool truthmatch = theAnalysis->McTruthMatch(antiproton[j]);
			ntpAntiProton->Column("MCTruthMatch", (bool) truthmatch);

			int gemhit = GemHits(antiproton[j]);

			int count = 0;
			if (moth==88888 && gemhit==1 && truthmatch==1) count=1;

			ntpAntiProton->Column("GemHit", (int) count, 0);


	    }

	    ntpAntiProton->DumpData();

	  }


	  //Write output
	  out->cd();

	  ntpPiMinus ->GetInternalTree()->Write();
	  ntpPiPlus->GetInternalTree()->Write();
	  ntpKaonMinus ->GetInternalTree()->Write();
	  ntpKaonPlus->GetInternalTree()->Write();
	  ntpProton->GetInternalTree()->Write();
	  ntpAntiProton->GetInternalTree()->Write();


	  out->Save();

	  timer.Stop();
	  Double_t rtime = timer.RealTime();
	  Double_t ctime = timer.CpuTime();

	  cout<<"Macro finisched successfully."<<endl;
	  cout<<"Realtime: "<<rtime<<" s, CPU time: "<<ctime<<" s"<<endl;
	  cout<<endl;


	  exit(0);

	}

//.........这里部分代码省略.........
	iGenJetB3=iJet;
	genJetB3 = (Jet*) branchGenJet->At(iGenJetB3);
      }
      else if ((jetB4) && (deltaR(genJet->Eta, jetB4->Eta, genJet->Phi, jetB4->Phi) < MAX_MATCH_DIST) ) {
	iGenJetB4=iJet;
	genJetB4 = (Jet*) branchGenJet->At(iGenJetB4);
      }
      else if ((jet1) && (deltaR(genJet->Eta, jet1->Eta, genJet->Phi, jet1->Phi) < MAX_MATCH_DIST) ) {
	iGenJetVBF1=iJet;
	genJetVBF1 = (Jet*) branchGenJet->At(iGenJetVBF1);
      }
      else if ((jet2) && (deltaR(genJet->Eta, jet2->Eta, genJet->Phi, jet2->Phi) < MAX_MATCH_DIST) ) {
	iGenJetVBF2=iJet;
	genJetVBF2 = (Jet*) branchGenJet->At(iGenJetVBF2);
      }
    }

    LorentzVector vGenJetB1(0,0,0,0);
    if (genJetB1) {
      vGenJetB1.SetPt(genJetB1->PT);
      vGenJetB1.SetEta(genJetB1->Eta);
      vGenJetB1.SetPhi(genJetB1->Phi);
      vGenJetB1.SetM(genJetB1->Mass);
      sGenJetB1 = &vGenJetB1;
    }
    else sGenJetB1 = &nothing;

    LorentzVector vGenJetB2(0,0,0,0);
    if (genJetB2) {
      vGenJetB2.SetPt(genJetB2->PT);
      vGenJetB2.SetEta(genJetB2->Eta);
      vGenJetB2.SetPhi(genJetB2->Phi);
      vGenJetB2.SetM(genJetB2->Mass);
      sGenJetB2 = &vGenJetB2;
    }
    else sGenJetB2 = &nothing;

    LorentzVector vGenJetB3(0,0,0,0);
    if (genJetB3) {
      vGenJetB3.SetPt(genJetB3->PT);
      vGenJetB3.SetEta(genJetB3->Eta);
      vGenJetB3.SetPhi(genJetB3->Phi);
      vGenJetB3.SetM(genJetB3->Mass);
      sGenJetB3 = &vGenJetB3;
    }
    else sGenJetB3 = &nothing;

    LorentzVector vGenJetB4(0,0,0,0);
    if (genJetB4) {
      vGenJetB4.SetPt(genJetB4->PT);
      vGenJetB4.SetEta(genJetB4->Eta);
      vGenJetB4.SetPhi(genJetB4->Phi);
      vGenJetB4.SetM(genJetB4->Mass);
      sGenJetB4 = &vGenJetB4;
    }
    else sGenJetB4 = &nothing;

    LorentzVector vGenJetVBF1(0,0,0,0);
    if (genJetVBF1) {
      vGenJetVBF1.SetPt(genJetVBF1->PT);
      vGenJetVBF1.SetEta(genJetVBF1->Eta);
      vGenJetVBF1.SetPhi(genJetVBF1->Phi);
      vGenJetVBF1.SetM(genJetVBF1->Mass);
      sGenJetVBF1 = &vGenJetVBF1;
    }
    else sGenJetVBF1 = &nothing;

    LorentzVector vGenJetVBF2(0,0,0,0);
    if (genJetVBF2) {
      vGenJetVBF2.SetPt(genJetVBF2->PT);
      vGenJetVBF2.SetEta(genJetVBF2->Eta);
      vGenJetVBF2.SetPhi(genJetVBF2->Phi);
      vGenJetVBF2.SetM(genJetVBF2->Mass);
      sGenJetVBF2 = &vGenJetVBF2;
    }
    else sGenJetVBF2 = &nothing;
    /*
    cout << "Gen Jets " << endl;
    if (genJetB1) cout << "1 " << genJetB1->PT << " " << genJetB1->Eta << endl;
    if (genJetB2) cout << "2 " << genJetB2->PT << " " << genJetB2->Eta << endl;
    if (genJetB3) cout << "3 " << genJetB3->PT << " " << genJetB3->Eta << endl;
    if (genJetB4) cout << "4 " << genJetB4->PT << " " << genJetB4->Eta << endl;
    if (genJetVBF1) cout << "V1 " << genJetVBF1->PT << " " << genJetVBF1->Eta << endl;
    if (genJetVBF2) cout << "V2 " << genJetVBF2->PT << " " << genJetVBF2->Eta << endl;
    cout << endl;
    */
    outTree->Fill();

  } // end event loop

  outFile->Write();
  outFile->Save();

  cout << endl;
  cout << "matched  : " << iMatched << endl;
  cout << "not      : " << iNot << endl;
  cout << "too much : " << iTwo << endl;
  cout << "total    : " << iMatched+iNot+iTwo << endl;

}

//.........这里部分代码省略.........
    tree->SetBranchStatus("m_b_DTF_jpsi_chic1_constr",1);
    tree->SetBranchStatus("dtf_b",1);
    tree->SetBranchStatus("c2ip_b",1);
    tree->SetBranchStatus("dls_b",1);
    tree->SetBranchStatus("c2ip_kminus",1);
    tree->SetBranchStatus("c2ip_pplus",1);
    tree->SetBranchStatus("minCl_gamma",1);
    tree->SetBranchStatus("m_chic",1);
    tree->SetBranchStatus("trig_b_l0tos_tos",1);
    tree->SetBranchStatus("trig_b_l1tos_tos",1);
    tree->SetBranchStatus("trig_b_l2tos_tos",1);
    tree->SetBranchStatus("trgh_track",1);
    tree->SetBranchStatus("ann_kaon",1);
    tree->SetBranchStatus("K_ProbNNk",1);
    tree->SetBranchStatus("ann_proton",1);
    tree->SetBranchStatus("P_ProbNNp",1);
    tree->SetBranchStatus("ann_mu",1);
    tree->SetBranchStatus("minann_mu",1);
    tree->SetBranchStatus("m_chicp",1);
    tree->SetBranchStatus("m_pK",1);
    tree->SetBranchStatus("m_jpsi",1);

    // -- this file is just here to make the 'CopyTree' happy
    TFile* dummyFile = new TFile("~/cern/ntuples/dummy.root","RECREATE");
    TTree* rTree1 = tree->CopyTree( cuts.c_str() );

    float mass, mass_chicp, mass_pK, mass_jpsi, mass_chic, dtf_b, c2ip_b, dls_b, c2ip_kminus, c2ip_pplus, minCl_gamma, trgh_track, ann_kaon, K_ProbNNk, ann_proton, P_ProbNNp, ann_mu, minann_mu;

    float trig_b_l0tos_tos, trig_b_l1tos_tos, trig_b_l2tos_tos;

    rTree1->SetBranchAddress("m_b_DTF_jpsi_chic1_constr", &mass);
    rTree1->SetBranchAddress("m_chicp", &mass_chicp);
    rTree1->SetBranchAddress("m_pK", &mass_pK);
    rTree1->SetBranchAddress("m_jpsi", &mass_jpsi);
    rTree1->SetBranchAddress("m_chic", &mass_chic);
    rTree1->SetBranchAddress("dtf_b", &dtf_b);
    rTree1->SetBranchAddress("c2ip_b", &c2ip_b);
    rTree1->SetBranchAddress("dls_b", &dls_b);
    rTree1->SetBranchAddress("c2ip_kminus", &c2ip_kminus);
    rTree1->SetBranchAddress("c2ip_pplus", &c2ip_pplus);
    rTree1->SetBranchAddress("minCl_gamma", &minCl_gamma);
    rTree1->SetBranchAddress("trig_b_l0tos_tos", &trig_b_l0tos_tos);
    rTree1->SetBranchAddress("trig_b_l1tos_tos", &trig_b_l1tos_tos);
    rTree1->SetBranchAddress("trig_b_l2tos_tos", &trig_b_l2tos_tos);
    rTree1->SetBranchAddress("trgh_track", &trgh_track);
    rTree1->SetBranchAddress("ann_kaon", &ann_kaon);
    rTree1->SetBranchAddress("K_ProbNNk", &K_ProbNNk);
    rTree1->SetBranchAddress("ann_proton", &ann_proton);
    rTree1->SetBranchAddress("P_ProbNNp", &P_ProbNNp);
    rTree1->SetBranchAddress("ann_mu", &ann_mu);
    rTree1->SetBranchAddress("minann_mu", &minann_mu);


    TFile* rFile = new TFile( outFilename.c_str() ,"RECREATE");
    TTree* rTree2 = new TTree();
    //rTree2->SetName("sideband");
    rTree2->SetName("signal");
    rTree2->Branch("mass", &mass, "mass");
    rTree2->Branch("mass_chicp", &mass_chicp, "mass_chicp/F");
    rTree2->Branch("mass_pK", &mass_pK, "mass_pK/F");
    rTree2->Branch("mass_jpsi", &mass_jpsi, "mass_jpsi/F");
    rTree2->Branch("mass_chic", &mass_chic, "mass_chic/F");
    rTree2->Branch("dtf_b", &dtf_b, "dtf_b/F");
    rTree2->Branch("c2ip_b", &c2ip_b, "c2ip_b/F");
    rTree2->Branch("dls_b", &dls_b, "dls_b/F");
    rTree2->Branch("c2ip_kminus", &c2ip_kminus, "c2ip_kminus/F");
    rTree2->Branch("c2ip_pplus", &c2ip_pplus, "c2ip_pplus/F");
    rTree2->Branch("minCl_gamma", &minCl_gamma, "minCl_gamma/F");
    rTree2->Branch("trig_b_l0tos_tos", &trig_b_l0tos_tos, "trig_b_l0tos_tos/F");
    rTree2->Branch("trig_b_l1tos_tos", &trig_b_l1tos_tos, "trig_b_l1tos_tos/F");
    rTree2->Branch("trig_b_l2tos_tos", &trig_b_l2tos_tos, "trig_b_l2tos_tos/F");
    rTree2->Branch("trgh_track", &trgh_track, "trgh_track/F");
    rTree2->Branch("ann_kaon", &ann_kaon, "ann_kaon/F");
    rTree2->Branch("K_ProbNNk", &K_ProbNNk, "K_ProbNNk/F");
    rTree2->Branch("ann_proton", &ann_proton, "ann_proton/F");
    rTree2->Branch("P_ProbNNp", &P_ProbNNp, "P_ProbNNp/F");
    rTree2->Branch("ann_mu", &ann_mu, "ann_mu/F");
    rTree2->Branch("minann_mu", &minann_mu, "minann_mu/F");


    int percentCounter = 1;

    for(int i = 0; i < rTree1->GetEntries(); ++i) {

        const int percent = (int)(rTree1->GetEntries()/100.0);

        if( i == percent*percentCounter ) {
            //std::cout << percentCounter << " %" << std::endl;
            percentCounter++;
        }

        rTree1->GetEntry(i);
        rTree2->Fill();
        //std::cout << mass << std::endl;
    }
    rTree2->Print();
    rTree2->Write();
    rFile->Save();

}

void effCalc(){ 

    // -- define tuple file name, tuple name and cuts to apply
    // -- and also the name of the output file
    
    const std::string filename = "/afs/cern.ch/work/a/apmorris/public/Lb2chicpK/signal_samples/Lb2chicpK_MC_2011_2012_signal_withbdt.root";
    //const std::string outFilename("/afs/cern.ch/work/a/apmorris/public/Lb2chicpK/signal_samples/");
        
    const std::string treename("withbdt");
    
    const std::string outfile = "~/cern/new_plots/signal_MC_gamma_PT_cut.pdf";
    const std::string xaxis   = "m(#chi_{c1}pK^{-}) (MeV/c^{2})";
    const std::string yaxis   = "Events / 2 MeV/c^{2}";
    const std::string cuts("gamma_PT > 500. ");
    
    //TStopwatch sw;
    //sw.Start();
    std::cout << "Opening file: " << filename.c_str() << endl;
    
    TFile* file = TFile::Open( filename.c_str() );
    if( !file ) std::cout << "file " << filename << " does not exist" << std::endl;

    TTree* tree = (TTree*)file->Get( treename.c_str() );
    if( !tree ) std::cout << "tree " << treename << " does not exist" << std::endl;
  
    // -- activate the branches you need
  
    tree->SetBranchStatus("*", 1);   
      
    //double n_pre = tree->GetEntries("bdtg>=0.85");
    //std::cout << "# of events in original tree = " << n_pre << endl;  
      
    // -- this file is just here to make the 'CopyTree' happy
    
    std::cout << "Copying tree: " << treename.c_str() << endl;
    
    TFile* newFile = new TFile("/afs/cern.ch/work/a/apmorris/public/Lb2chicpK/signal_samples/dummy.root","RECREATE");
    TTree* rTree1 = tree->CopyTree( cuts.c_str() );
    
    //rTree1->Print();
    rTree1->Write();
    newFile->Save();

    
    double n_post = rTree1->GetEntries();
    
    // Including the LHCb name, formatting
    lhcbName = new TPaveText(gStyle->GetPadLeftMargin() + 0.65,
            0.87 - gStyle->GetPadTopMargin(),
            gStyle->GetPadLeftMargin() + 0.75,
            0.95 - gStyle->GetPadTopMargin(),
            "BRNDC");
    lhcbName->AddText("LHCb");
    lhcbName->SetFillColor(0);
    lhcbName->SetTextAlign(12);
    lhcbName->SetBorderSize(0);

    gStyle->SetOptStat(0);	

/*
    titleBox = new TPaveText(gStyle->GetPadLeftMargin() + 0.3,
            0.98 - gStyle->GetPadTopMargin(),
            0.98 - gStyle->GetPadLeftMargin(),
            0.83 - gStyle->GetPadTopMargin(),
            "BRNDC");
    titleBox->AddText(title.c_str());
    titleBox->SetFillColor(0);
    titleBox->SetTextAlign(12);
    titleBox->SetBorderSize(0);
*/
    
    

    // plotting ------------------------------------------------------------------- 

    // Define the histograms
    TH1D* h1 = new TH1D ("h1", "Effect of cutting #gamma pT on simulated sample post-selection", 100, 5500., 5700.);
    TH1D* h2 = new TH1D ("h2", "Effect of cutting #gamma pT on simulated sample post-selection", 100, 5500., 5700.);
    
    
    // Differentiating the two histos from each other
    //h1->Sumw2(); 	
    //h2->Sumw2();
    
    //h1->SetMarkerStyle(kFullDotLarge);
    h1->SetMarkerSize(0.7);
   
    h2->SetMarkerColor(kRed);
    //h2->SetMarkerStyle(kFullDotLarge);
    h2->SetMarkerSize(0.7);
    h2->SetLineColor(kRed);
    

    // Draw the histograms 
    rTree1->Draw("Lambda_b0_DTF_MASS_constr1>>h1", "", "goff"); 	
    tree->Draw("Lambda_b0_DTF_MASS_constr1>>h2", "", "goff");

    
    // Setting up canvas
    TCanvas* c1 = new TCanvas();
//.........这里部分代码省略.........

void computeSummary(string fileName_="",string payloadName_=""){
   if (fileName_==""){cout<<"No file specified."<<endl; return;}
   TFile * f = new TFile(saveFile.c_str(),"UPDATE");
   TTree * t = (TTree*)f->Get("payloadSummary"); 
   if (!t)return;
   float  MPV[5];
   float  Gain[5];
   gROOT->ProcessLine("#include <iostream>");
   string *payloadName = new string;
   payloadName->assign(payloadName_); 
   t->SetBranchAddress("MPV"        ,&MPV[0]           );
   t->SetBranchAddress("Gain"       ,&Gain[0]          );
   t->SetBranchAddress("MPV_TIB"    ,&MPV[1]           );
   t->SetBranchAddress("MPV_TID"    ,&MPV[2]           );
   t->SetBranchAddress("MPV_TOB"    ,&MPV[3]           );
   t->SetBranchAddress("MPV_TEC"    ,&MPV[4]           );
   t->SetBranchAddress("Gain_TIB"   ,&Gain[1]          );
   t->SetBranchAddress("Gain_TID"   ,&Gain[2]          );
   t->SetBranchAddress("Gain_TOB"   ,&Gain[3]          );
   t->SetBranchAddress("Gain_TEC"   ,&Gain[4]          );
   t->SetBranchAddress("PayloadName",&payloadName      );
   
   float NAPV[5];
   for(int i=0; i<5; i++){MPV[i]=0; Gain[i]=0; NAPV[i]=0;}
   
   TFile * f2 = new TFile(fileName_.c_str());
   TTree * t2 = (TTree*) f2->Get("SiStripCalib/APVGain");
   float MPV_origin; double Gain_origin; UChar_t subDet; bool isMasked;
   t2->SetBranchAddress("FitMPV"  ,&MPV_origin );
   t2->SetBranchAddress("Gain"    ,&Gain_origin);
   t2->SetBranchAddress("SubDet"  ,&subDet     );
   t2->SetBranchAddress("isMasked",&isMasked   );

   //Looping over tree...
   printf("Progressing Bar              :0%%       20%%       40%%       60%%       80%%       100%%\n");
   printf("Looping on the Tree          :");
   int TreeStep = t2->GetEntries()/50;if(TreeStep==0)TreeStep=1;
   for (unsigned int ientry = 0; ientry < t2->GetEntries(); ientry++) {
      if(ientry%TreeStep==0){printf(".");fflush(stdout);}
      t2->GetEntry(ientry);

      if(subDet>2 && Gain_origin > 0 && MPV_origin > 0 && ! isMasked){
         NAPV[subDet-2]++;
         NAPV[0]++;
         Gain[0]+=Gain_origin;
         MPV[0]+=MPV_origin;
         Gain[subDet-2]+=Gain_origin;
         MPV [subDet-2]+=MPV_origin;
      }

   }
   cout<<endl;
   f2->Close();
   
   for (int i=0; i<5; i++){MPV[i]/=NAPV[i]; Gain[i]/=NAPV[i];}

   f->cd();
   
   t->Fill();
   t->Write();
   f->Save();
   delete payloadName;
   f->Close();
}

//.........这里部分代码省略.........
		// *** and again with PidAlgoMvd;PidAlgoStt;PidAlgoDrc and loose selection
		theAnalysis->FillList(muplus,  "MuonLoosePlus",  "PidAlgoMvd;PidAlgoStt;PidAlgoDrc;PidAlgoMdtHardCuts");
		theAnalysis->FillList(muminus, "MuonLooseMinus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc;PidAlgoMdtHardCuts");
		theAnalysis->FillList(piplus,  "PionLoosePlus",  "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
		theAnalysis->FillList(piminus, "PionLooseMinus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
		
		jpsi.Combine(muplus, muminus);
		for (j=0;j<jpsi.GetLength();++j) hjpsim_lpid->Fill( jpsi[j]->M() );
		jpsi.Select(jpsiMassSel);
		
		psi2s.Combine(jpsi, piplus, piminus);
		for (j=0;j<psi2s.GetLength();++j) hpsim_lpid->Fill( psi2s[j]->M() );
		
		
		// ***
		// *** TIGHT PID combinatorics
		// ***
		
		// *** and again with PidAlgoMvd;PidAlgoStt and tight selection
		theAnalysis->FillList(muplus,  "MuonTightPlus",  "PidAlgoMdtHardCuts");
		theAnalysis->FillList(muminus, "MuonTightMinus", "PidAlgoMdtHardCuts");
		theAnalysis->FillList(piplus,  "PionLoosePlus",  "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
		theAnalysis->FillList(piminus, "PionLooseMinus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
		
		jpsi.Combine(muplus, muminus);
		for (j=0;j<jpsi.GetLength();++j) hjpsim_tpid->Fill( jpsi[j]->M() );
		jpsi.Select(jpsiMassSel);
		
		psi2s.Combine(jpsi, piplus, piminus);
		for (j=0;j<psi2s.GetLength();++j) hpsim_tpid->Fill( psi2s[j]->M() );
		
	}
		
	// *** write out all the histos
	out->cd();
	
	hmomtrk->Write();
	hthttrk->Write();
	
	hjpsim_all->Write();
	hpsim_all->Write();
	hjpsim_lpid->Write();
	hpsim_lpid->Write();
	hjpsim_tpid->Write();
	hpsim_tpid->Write();
	hjpsim_trpid->Write();
	hpsim_trpid->Write();
	
	hjpsim_ftm->Write();
	hpsim_ftm->Write();
	hjpsim_nm->Write();
	hpsim_nm->Write();
	
	hpsim_diff->Write();
	hjpsim_diff->Write();
	
	hjpsim_vf->Write();
	hjpsim_4cf->Write();
	hjpsim_mcf->Write();
	
	hjpsi_chi2_vf->Write();
	hpsi_chi2_4c->Write();
	hjpsi_chi2_mf->Write();
			
	hjpsi_prob_vf->Write();
	hpsi_prob_4c->Write();
	hjpsi_prob_mf->Write();
			
	hvpos->Write();
		
	out->Save();
        
	// Extract the maximal used memory an add is as Dart measurement
	// This line is filtered by CTest and the value send to CDash
	FairSystemInfo sysInfo;
	Float_t maxMemory=sysInfo.GetMaxMemory();
	cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">";
	cout << maxMemory;
	cout << "</DartMeasurement>" << endl;
  
	fTimer.Stop();
	Double_t rtime = fTimer.RealTime();
	Double_t ctime = fTimer.CpuTime();
  
	Float_t cpuUsage=ctime/rtime;
	cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
	cout << cpuUsage;
	cout << "</DartMeasurement>" << endl;
  
	cout << endl;
	cout << "Real time " << rtime << " s, CPU time " << ctime
		  << "s" << endl;
	cout << "CPU usage " << cpuUsage*100. << "%" << endl;
	cout << "Max Memory " << maxMemory << " MB" << endl;
   
	cout << "Macro finished successfully." << endl;

        exit(0);
	
}

//.........这里部分代码省略.........
            z2 = (TGenParticle*) ((*part)[iZ2]);

            z1_pt=z1->pt;
            z1_eta=z1->eta;
            z1_phi=z1->phi;
            z1_mass=z1->mass;

            z2_pt=z2->pt;
            z2_eta=z2->eta;
            z2_phi=z2->phi;
            z2_mass=z2->mass;
        }

        if(iZ!=-1){
            const TGenParticle *parZ = (TGenParticle*) ((*part)[iZ]);

            z_pt = parZ->pt;
            z_eta = parZ->eta;
            z_phi = parZ->phi;
            z_mass = parZ->mass;
        }

        recoz_pt = vRecoZ.Pt();
        recoz_eta = vRecoZ.Eta();
        recoz_phi = vRecoZ.Phi();
        recoz_mass = vRecoZ.M();

        if(iT1!=-1 && iT2!=-1){
            const TGenJet *taujet1, *taujet2;

            taujet1 = (TGenJet*) ((*jet)[iT1]);
            taujet2 = (TGenJet*) ((*jet)[iT2]);

            jetTau1_pt=taujet1->pt;
            jetTau1_eta=taujet1->eta;
            jetTau1_phi=taujet1->phi;
            jetTau1_mass=taujet1->mass;

            jetTau2_pt=taujet2->pt;
            jetTau2_eta=taujet2->eta;
            jetTau2_phi=taujet2->phi;
            jetTau2_mass=taujet2->mass;
        } 

        genTau1_pt=genTau1->pt;
        genTau1_eta=genTau1->eta;
        genTau1_phi=genTau1->phi;
        genTau1_mass=genTau1->mass;

        genTau2_pt=genTau2->pt;
        genTau2_eta=genTau2->eta;
        genTau2_phi=genTau2->phi;
        genTau2_mass=genTau2->mass;

        visTau1_pt=vVisTau1.Pt();
        visTau1_eta=vVisTau1.Eta();
        visTau1_phi=vVisTau1.Phi();
        visTau1_mass=vVisTau1.M();

        visTau2_pt=vVisTau2.Pt();
        visTau2_eta=vVisTau2.Eta();
        visTau2_phi=vVisTau2.Phi();
        visTau2_mass=vVisTau2.M();

        cpions1_pt=cpions1.Pt();
        cpions1_eta=cpions1.Eta();
        cpions1_phi=cpions1.Phi();
        cpions1_mass=cpions1.M();

        cpions2_pt=cpions2.Pt();
        cpions2_eta=cpions2.Eta();
        cpions2_phi=cpions2.Phi();
        cpions2_mass=cpions2.M();

        npions1_pt=npions1.Pt();
        npions1_eta=npions1.Eta();
        npions1_phi=npions1.Phi();
        npions1_mass=npions1.M();

        npions2_pt=npions2.Pt();
        npions2_eta=npions2.Eta();
        npions2_phi=npions2.Phi();
        npions2_mass=npions2.M();

        for (Int_t j=0; j<part->GetEntries(); j++)
        {
            const TGenParticle* genloop = (TGenParticle*) ((*part)[j]);
            if ( genloop->pt<20) continue;
            if ( abs(genloop->pdgId)!=13 && abs(genloop->pdgId)!=11 ) continue;
            nLeptons++;
        }

        outtree->Fill();

    }

    outfile->Write();
    outfile->Save();

}

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

    //hhMVA::MVAType why1=hhMVA::kTauTau;
    hhMVA::MVAType why2=hhMVA::kMuTau;
    hhMVA::MVAType why3=hhMVA::kElTau;
    hhMVA::MVAType why4=hhMVA::kElMu;

    //hhMVA *ttMVA = new hhMVA();
    hhMVA *mtMVA = new hhMVA();
    hhMVA *etMVA = new hhMVA();
    hhMVA *emMVA = new hhMVA();

    //ttMVA->Intialize(why1);
    mtMVA->Intialize(why2);
    etMVA->Intialize(why3);
    emMVA->Intialize(why4);

    for (Int_t i=0; i<eventChain->GetEntries(); i++) {
        eventChain->GetEntry(i);

        eventWeight/=float(nevents);
        bdtVal=999;

        // hack to fix modified SM samples
        //if (sampleNo>95 && sampleNo<100) eventWeight=2.92/float(nevents);

        if (isBBTT!=1) continue;
        if (ptB1<30 || ptB2<30) continue;
        //if (TMath::Sqrt((etaTau1-etaTau2)*(etaTau1-etaTau2)+(phiTau1-phiTau2)*(phiTau1-phiTau2))<0.4) continue;
        /*
        TLorentzVector v1;
        TLorentzVector v2;

        v1.SetPtEtaPhiM(ptB1, etaB1, phiB1, mB1);
        v2.SetPtEtaPhiM(ptB2, etaB2, phiB2, mB2);

        TMatrixD m1(3,3);
        TMatrixD m2(3,3);

        m1.Zero();
        m2.Zero();
        m1(0,0) = ErrEt (v1.Et(), v1.Eta()); // et
        m1(1,1) = ErrEta(v1.Et(), v1.Eta()); // eta
        m1(2,2) = ErrPhi(v1.Et(), v1.Eta()); // phi
        m2(0,0) = ErrEt (v2.Et(), v2.Eta()); // et
        m2(1,1) = ErrEta(v2.Et(), v2.Eta()); // eta
        m2(2,2) = ErrPhi(v2.Et(), v2.Eta()); // phi

        TFitParticleEtEtaPhi *jet1 = new TFitParticleEtEtaPhi( "Jet1", "Jet1", &v1, &m1 );
        TFitParticleEtEtaPhi *jet2 = new TFitParticleEtEtaPhi( "Jet2", "Jet2", &v2, &m2 );

        TFitConstraintM *mCons1 = new TFitConstraintM( "HMassConstraint", "HMass-Constraint", 0, 0 , 125.);
        mCons1->addParticles1( jet1, jet2 );

        TKinFitter* fitter = new TKinFitter("fitter", "fitter");
        fitter->addMeasParticle( jet1 );
        fitter->addMeasParticle( jet2 );
        fitter->addConstraint( mCons1 );

        //Set convergence criteria
        fitter->setMaxNbIter( 30 );
        fitter->setMaxDeltaS( 1e-2 );
        fitter->setMaxF( 1e-1 );
        fitter->setVerbosity(1);

        TVector2 b1_i; b1_i.SetMagPhi(fitter->get4Vec(0)->Pt(), fitter->get4Vec(0)->Phi());
        TVector2 b2_i; b2_i.SetMagPhi(fitter->get4Vec(1)->Pt(), fitter->get4Vec(1)->Phi());

        //Perform the fit
        conv=fitter->fit();

        chi2=fitter->getS();

        TVector2 metTemp; metTemp.SetMagPhi(pileupMet,pileupMetPhi);
        TVector2 b1_f; b1_f.SetMagPhi(fitter->get4Vec(0)->Pt(), fitter->get4Vec(0)->Phi());
        TVector2 b2_f; b2_f.SetMagPhi(fitter->get4Vec(1)->Pt(), fitter->get4Vec(1)->Phi());

        metTemp=metTemp+b1_i+b2_i-b1_f-b2_f;
        corrMet=metTemp.Mod();
        corrMetPhi=metTemp.Phi();
        */
        if (tauCat1==1 && tauCat2==1 && ptTau1>45 && ptTau2>45) {
            bdtVal=999;//ttMVA->GetBDTValue(mTT, ptTT, mBB1, ptBB1, dRBB1, dRTT);
        }
        else if (tauCat1==1 && tauCat2==3 && ptTau1>30 && ptTau2>20) {
            bdtVal=mtMVA->GetBDTValue(mTT, ptTT, mBB1, ptBB1, mHH, ptHH, mt2pileup, dRBB1, dRTT, dRHH);
        }
        else if (tauCat1==1 && tauCat2==2 && ptTau1>30 && ptTau2>20) {
            bdtVal=etMVA->GetBDTValue(mTT, ptTT, mBB1, ptBB1, mHH, ptHH, mt2pileup, dRBB1, dRTT, dRHH);
        }
        else if (tauCat1==3 && tauCat2==2 && ptTau1>20 && ptTau2>20) {
            bdtVal=emMVA->GetBDTValue(mTT, ptTT, mBB1, ptBB1, mHH, ptHH, mt2pileup, dRBB1, dRTT, dRHH);
        }
        else continue;
        outTree->Fill();
    }

    outFile->Write();
    outFile->Save();

}

//.........这里部分代码省略.........
  if (!inpF || !inpF->IsOpen()) {
    std::cerr << "[ERROR]: Could not open input file: " << InpFileName
              << std::endl;
    exit(1);
  }

  TTree *stdhep = static_cast<TTree *>(inpF->Get("giRooTracker"));

  if (!stdhep) {
    std::cerr
        << "[ERROR]: Could not read TTree (\"giRooTracker\") from input file: "
        << InpFileName << std::endl;
    exit(1);
  }

  Double_t EvtWght;
  Int_t GiStdHepN;
  Int_t GiStdHepPdg[100];
  Int_t GiStdHepStatus[100];
  Double_t GiStdHepP4[100][4];

  stdhep->SetBranchAddress("EvtWght", &EvtWght);
  stdhep->SetBranchAddress("StdHepN", &GiStdHepN);
  stdhep->SetBranchAddress("StdHepPdg", GiStdHepPdg);
  stdhep->SetBranchAddress("StdHepP4", GiStdHepP4);
  stdhep->SetBranchAddress("StdHepStatus", GiStdHepStatus);

  TFile *oupF = new TFile(OupFileName, "RECREATE");
  TTree *oupT = new TTree("CC1PipQ2", "");
  Float_t Q2;
  oupT->Branch("Q2", &Q2);
  oupT->Branch("EvtWght", &EvtWght);

  Long64_t nentries = stdhep->GetEntries();
  for (Long64_t evt = 0; evt < nentries; ++evt) {
    stdhep->GetEntry(evt);

    TLorentzVector pnu(0, 0, 0, 0);
    TLorentzVector pmu(0, 0, 0, 0);

    Int_t NPiPlus = 0;
    Int_t NOtherPi = 0;

    if (db) {
      std::cout << "Ev[" << evt << "] ------- " << std::endl;
    }

    // Loop through particle stack
    for (Int_t prt = 0; prt < GiStdHepN; ++prt) {
      // Inital numu
      if ((GiStdHepStatus[prt] == 0) && (GiStdHepPdg[prt] == 14)) {
        pnu = TLorentzVector(GiStdHepP4[prt][0], GiStdHepP4[prt][1],
                             GiStdHepP4[prt][2], GiStdHepP4[prt][3]);
        if (db) {
          std::cout << "\tFound nu at " << prt << " Mom: ("
                    << GiStdHepP4[prt][0] << ", " << GiStdHepP4[prt][1] << ", "
                    << GiStdHepP4[prt][2] << ", " << GiStdHepP4[prt][3] << ") "
                    << std::endl;
        }
      }
      // Final mu
      if ((GiStdHepStatus[prt] == 1) && (GiStdHepPdg[prt] == 13)) {
        pmu = TLorentzVector(GiStdHepP4[prt][0], GiStdHepP4[prt][1],
                             GiStdHepP4[prt][2], GiStdHepP4[prt][3]);
        if (db) {
          std::cout << "\tFound mu at " << prt << " Mom: ("
                    << GiStdHepP4[prt][0] << ", " << GiStdHepP4[prt][1] << ", "
                    << GiStdHepP4[prt][2] << ", " << GiStdHepP4[prt][3] << ") "
                    << std::endl;
        }
      }
      // Final pi+
      if ((GiStdHepStatus[prt] == 1) && (GiStdHepPdg[prt] == 211)) {
        NPiPlus++;
      }
      // Final other pi
      if ((GiStdHepStatus[prt] == 1) &&
          ((GiStdHepPdg[prt] == 111) || (GiStdHepPdg[prt] == -211))) {
        NOtherPi++;
      }
    }
    // cc1pip selection
    if ((NPiPlus == 1) && (NOtherPi == 0) && (pnu.Vect().Mag2() > 0) &&
        (pmu.Vect().Mag2() > 0)) {
      Q2 = -1. * (pnu - pmu).Mag2();
      if (db) {
        std::cout << "Ev[" << evt << "] -- Q2: " << Q2
                  << ", pnu.Mag(): " << pnu.Vect().Mag()
                  << ", pmu.Mag(): " << pmu.Vect().Mag() << std::endl;
      }
      oupT->Fill();
    }
    if (db) {
      std::cout << "====================" << std::endl;
    }
  }
  oupT->Write();
  oupF->Write();
  oupF->Save();
}

//.........这里部分代码省略.........
    rTree2->Branch("Lambda_b0_L0DiMuonDecision_TOS",&Lambda_b0_L0DiMuonDecision_TOS,"Lambda_b0_L0DiMuonDecision_TOS/B");
    rTree2->Branch("Lambda_b0_L0MuonDecision_TOS",  &Lambda_b0_L0MuonDecision_TOS,  "Lambda_b0_L0MuonDecision_TOS/B");
    rTree2->Branch("Lambda_b0_FDS",                 &Lambda_b0_FDS,                 "Lambda_b0_FDS/D");  
    
    rTree2->Branch("Lambda_b0_Hlt1DiMuonHighMassDecision_TOS",      &Lambda_b0_Hlt1DiMuonHighMassDecision_TOS,      "Lambda_b0_Hlt1DiMuonHighMassDecision_TOS/B");  
    rTree2->Branch("Lambda_b0_Hlt1DiMuonLowMassDecision_TOS",       &Lambda_b0_Hlt1DiMuonLowMassDecision_TOS,       "Lambda_b0_Hlt1DiMuonLowMassDecision_TOS/B");  
    rTree2->Branch("Lambda_b0_Hlt1TrackMuonDecision_TOS",           &Lambda_b0_Hlt1TrackMuonDecision_TOS,           "Lambda_b0_Hlt1TrackMuonDecision_TOS/B");  
    rTree2->Branch("Lambda_b0_Hlt1TrackAllL0Decision_TOS",          &Lambda_b0_Hlt1TrackAllL0Decision_TOS,          "Lambda_b0_Hlt1TrackAllL0Decision_TOS/B");
    rTree2->Branch("Lambda_b0_Hlt1SingleMuonHighPTDecision_TOS",    &Lambda_b0_Hlt1SingleMuonHighPTDecision_TOS,    "Lambda_b0_Hlt1SingleMuonHighPTDecision_TOS/B");  
    rTree2->Branch("Lambda_b0_Hlt2DiMuonDetachedDecision_TOS",      &Lambda_b0_Hlt2DiMuonDetachedDecision_TOS,      "Lambda_b0_Hlt2DiMuonDetachedDecision_TOS/B"); 
    rTree2->Branch("Lambda_b0_Hlt2DiMuonDetachedJPsiDecision_TOS",  &Lambda_b0_Hlt2DiMuonDetachedJPsiDecision_TOS,  "Lambda_b0_Hlt2DiMuonDetachedJPsiDecision_TOS/B");
    rTree2->Branch("Lambda_b0_Hlt2DiMuonDetachedHeavyDecision_TOS", &Lambda_b0_Hlt2DiMuonDetachedHeavyDecision_TOS, "Lambda_b0_Hlt2DiMuonDetachedHeavyDecision_TOS/B");
    
    rTree2->Branch("Lambda_b0_pi0veto", &Lambda_b0_pi0veto, "Lambda_b0_pi0veto/D");
    rTree2->Branch("Lambda_b0_PT",      &Lambda_b0_PT,      "Lambda_b0_PT/D");
                            
                                              
    //for MC only                                                                   //
    rTree2->Branch("chi_c_BKGCAT",      &chi_c_BKGCAT,      "chi_c_BKGCAT/I");      //
    rTree2->Branch("Jpsi_BKGCAT",       &Jpsi_BKGCAT,       "Jpsi_BKGCAT/I");       //
    rTree2->Branch("Lambda_b0_BKGCAT",  &Lambda_b0_BKGCAT,  "Lambda_b0_BKGCAT/I");  //
    
    
//-----------------------------------------------------------------------------------    
    
	const double mK(493.677);
    const double mpi(139.57);
    const double mp(938.27);
    
    double m_pK, m_chicp;
    double proton_as_kaon_M,        proton_as_pion_M,           proton_as_proton_M;
    double m_chicpK_proton_as_kaon, m_chicpK_proton_as_pion,    m_chicpK_proton_as_proton; 
    double m_JpsipK;
    
    rTree2->Branch("m_JpsipK",  &m_JpsipK,  "m_JpsipK/D");
	rTree2->Branch("m_chicp",   &m_chicp,   "m_chicp/D");
    rTree2->Branch("m_pK",      &m_pK,      "m_pK/D");
    
    rTree2->Branch("proton_as_kaon_M",      &proton_as_kaon_M,      "proton_as_kaon_M/D");
    rTree2->Branch("proton_as_pion_M",      &proton_as_pion_M,      "proton_as_pion_M/D");
    rTree2->Branch("proton_as_proton_M",    &proton_as_proton_M,    "proton_as_proton_M/D");
    
    rTree2->Branch("m_chicpK_proton_as_kaon",     &m_chicpK_proton_as_kaon,     "m_chicpK_proton_as_kaon/D");
    rTree2->Branch("m_chicpK_proton_as_pion",     &m_chicpK_proton_as_pion,     "m_chicpK_proton_as_pion/D");
    rTree2->Branch("m_chicpK_proton_as_proton",   &m_chicpK_proton_as_proton,   "m_chicpK_proton_as_proton/D");
	
    int percentCounter = 1;

    for(int i = 0; i < rTree1->GetEntries(); ++i){
  
      const int percent = (int)(rTree1->GetEntries()/100.0);
    
      if( i == percent*percentCounter ){
        std::cout << percentCounter << " %" << std::endl;
        percentCounter++;
      }
      
      rTree1->GetEntry(i);
      
      double proton_P = sqrt(proton_PX*proton_PX + proton_PY*proton_PY + proton_PZ*proton_PZ) ;
      
      TLorentzVector proton_as_kaon(    proton_PX,  proton_PY,  proton_PZ,  sqrt(proton_P*proton_P + mK*mK));
      TLorentzVector proton_as_pion(    proton_PX,  proton_PY,  proton_PZ,  sqrt(proton_P*proton_P + mpi*mpi));
      TLorentzVector proton_as_proton(  proton_PX,  proton_PY,  proton_PZ,  sqrt(proton_P*proton_P + proton_M*proton_M));
      
      TLorentzVector kaon(  kaon_PX,    kaon_PY,    kaon_PZ,    kaon_PE);
      TLorentzVector chic(  chi_c_PX,   chi_c_PY,   chi_c_PZ,   chi_c_PE);
      TLorentzVector proton(proton_PX,  proton_PY,  proton_PZ,  proton_PE);
      TLorentzVector Jpsi(  Jpsi_PX,    Jpsi_PY,    Jpsi_PZ,    Jpsi_PE);
      
      proton_as_kaon_M            = proton_as_kaon.M();
      proton_as_pion_M            = proton_as_pion.M();
      proton_as_proton_M          = proton_as_proton.M();
        
      TLorentzVector chic_PasK_K  = chic + kaon + proton_as_kaon;
      TLorentzVector chic_PasPi_K = chic + kaon + proton_as_pion;
      TLorentzVector chic_PasP_K  = chic + kaon + proton_as_proton;
        
      m_chicpK_proton_as_kaon     = chic_PasK_K.M();
      m_chicpK_proton_as_pion     = chic_PasPi_K.M();
      m_chicpK_proton_as_proton   = chic_PasP_K.M();
      
      TLorentzVector JpsipK       = Jpsi   + proton + kaon;
      TLorentzVector chicp        = chic   + proton;
      TLorentzVector pK           = proton + kaon;
      
      m_JpsipK  = JpsipK.M();  
      m_chicp   = chicp.M();
      m_pK      = pK.M();
        
      rTree2->Fill();
         
    
    }
    
    rTree2->Print();
    rTree2->Write();
    rFile->Save();

}

//.........这里部分代码省略.........
    atuple[15]=k_disc_e;
    atuple[16]=k_drc_e;
    atuple[17]=k_mvd_e;
    atuple[18]=k_comb_e;

    // elec2
    PndPidProbability *drc_posi = (PndPidProbability *)drc_array->At(iy);
    PndPidProbability *disc_posi = (PndPidProbability *)disc_array->At(iy);
    PndPidProbability *mvd_posi = (PndPidProbability *)mvd_array->At(iy);
    PndPidProbability *stt_posi = (PndPidProbability *)stt_array->At(iy);
    PndPidProbability *emcb_posi = (PndPidProbability *)emcb_array->At(iy);
    Double_t k_drc_p = drc_posi->GetElectronPidProb();
    Double_t k_disc_p = disc_posi->GetElectronPidProb();
    Double_t k_mvd_p = mvd_posi->GetElectronPidProb();
    Double_t k_stt_p = stt_posi->GetElectronPidProb();
    Double_t k_emcb_p = emcb_posi->GetElectronPidProb();
    Double_t xx_p = (k_drc_p/(1-k_drc_p))*(k_disc_p/(1-k_disc_p))
      *(k_mvd_p/(1-k_mvd_p))*(k_stt_p/(1-k_stt_p))
      *(k_emcb_p/(1-k_emcb_p));
    Double_t k_comb_p = xx_p/(xx_p+1);

    atuple[19]=reTrack[1].P();
    atuple[20]=reTrack[1].Theta();
    atuple[21]=reTrack[1].Phi();
    atuple[22]=pc2->GetEmcRawEnergy();
    atuple[23]=pc2->GetEmcNumberOfCrystals();
    atuple[24]=k_emcb_p;
    atuple[25]=k_stt_p;
    atuple[26]=k_disc_p;
    atuple[27]=k_drc_p;
    atuple[28]=k_mvd_p;
    atuple[29]=k_comb_p;
    atuple[30]=bestTH;
    atuple[31]=bestPH;
    atuple[32]=bestCOST;

    NTev->Fill(atuple);
    if(j<5) {
      cout << "k_comb_e:" << k_comb_e<< " k_comb_p:" << k_comb_p << endl;
    }

    // standard cuts for histos
    Bool_t com_ele_1 = drc_ele->GetElectronPidProb() > 0.05 &&
      disc_ele->GetElectronPidProb() > 0.05 &&
      mvd_ele->GetElectronPidProb() > 0.05 &&
      stt_ele->GetElectronPidProb() > 0.05 &&
      emcb_ele->GetElectronPidProb() > 0.05;
    Bool_t com_ele_2 = k_comb_e > 0.9;
    Bool_t com_ele_3 = pc1->GetEmcNumberOfCrystals() > 5;
    Bool_t com_ele_4 = bestTH >=178. && bestTH <= 182.;
    Bool_t com_ele_5 = bestPH >=178. && bestPH <= 182.;
    Bool_t com_posi_1 = drc_posi->GetElectronPidProb() > 0.05 &&
      disc_posi->GetElectronPidProb() > 0.05 &&
      mvd_posi->GetElectronPidProb() > 0.05 &&
      stt_posi->GetElectronPidProb() > 0.05 &&
      emcb_posi->GetElectronPidProb() > 0.05;
    Bool_t com_posi_2 =  k_comb_p > 0.9;
    Bool_t com_posi_3 = pc2->GetEmcNumberOfCrystals() > 5;

    if(j<5) {
      cout << com_ele_1 << com_ele_2 << com_ele_3 << com_ele_4
	   << com_ele_4 << com_posi_1 << com_posi_2 << com_posi_3 << endl;
    }
    if (com_ele_1&&com_ele_2&&com_ele_3&&com_ele_4
        &&com_ele_5&&com_posi_1&&com_posi_2&&com_posi_3) {
      h_costheta_sel->Fill(bestCOST);

    }
  }  // loop j over events

  cout << " NEVcount: " << NEVcount << endl;
  h_efficiency->Divide(h_costheta_sel,h_costheta_mc,1,1,"B");

  // output
  cMA->cd(1);gPad->SetLogy(); h_costheta_mc->Draw();
  cMA->cd(2);gPad->SetLogy(); h_costheta->Draw();
  cMA->cd(3);gPad->SetLogy(); h_costheta_sel->Draw();
  cMA->cd(4);gPad->SetLogy(); h_efficiency->Draw();
  cMA->cd(5);gPad->SetLogy(); h_dtheta->Draw();
  cMA->cd(6);gPad->SetLogy(); h_dphi->Draw();
  cMA->cd(7); hpair->Draw();
  cMA->cd(0);

  out->cd();

  hpair->Write();
  h_costheta_mc->Write();
  h_costheta->Write();
  h_costheta_sel->Write();
  h_efficiency->Write();
  h_dtheta->Write();
  h_dphi->Write();

  NTev->Write();

  out->Save();

  cout << " Yahoo! " << endl;

}