C++ TChain::SetBranchAddress方法代码示例

//-------------------------------------------------------------------
void TestReadStuple()
{
	  TChain *myTree = new TChain("Stuple");

 //DATA
 //myTree->Add("/mnt/autofs/misc/nbay03.a/samantha/RESULTS/STUPLES/StupleV2_noMinVtx_noMinJet_TLphoEleRemoved.root");
 myTree->Add("/mnt/autofs/misc/nbay03.a/samantha/RESULTS/STUPLES/Stuple_noMinVtx_noMinJet_TLphoEleRemoved.root");

	Stuple *myStuple = new Stuple; 

	myTree->SetBranchStatus("pho_num", 1);		// number of electrons so we know how much to read from arrays
	myTree->SetBranchStatus("ele_num", 1);
	myTree->SetBranchStatus("jet_num", 1);

	myTree->SetBranchStatus("ele_Index", 1);
	myTree->SetBranchStatus("ele_Etc", 1);
	myTree->SetBranchStatus("ele_Ntight", 1);
	myTree->SetBranchStatus("ele_Nloose", 1);
	myTree->SetBranchStatus("ele_matchJetIndex", 1);	//if a matching jet is found and removed from jet list


	myTree->SetBranchAddress("pho_num", &myStuple->pho_num);		// number of electrons so we know how much to read from arrays
	myTree->SetBranchAddress("ele_num", &myStuple->ele_num);
	myTree->SetBranchAddress("jet_num", &myStuple->jet_num);

	myTree->SetBranchAddress("ele_Ntight", &myStuple->ele_Ntight);
	myTree->SetBranchAddress("ele_Nloose", &myStuple->ele_Nloose);
	myTree->SetBranchAddress("ele_Index", &myStuple->ele_Index);
	myTree->SetBranchAddress("ele_Etc", &myStuple->ele_Etc);
	myTree->SetBranchAddress("ele_matchJetIndex", &myStuple->ele_matchJetIndex);	//if a matching jet is found and removed from jet list


	//for (unsigned i =0 ; i < myTree->GetEntries(); ++i) {
	for (unsigned i =0 ; i < 100; ++i) {
		myTree->GetEntry(i);
		std::cout << myStuple->ele_num << std::endl;
		for (unsigned j = 0 ; j < myStuple->ele_num; ++j) {
			std::cout << myStuple->ele_matchJetIndex[j] << std::endl;

		}
	}
}

//.........这里部分代码省略.........
    std::cout << " -- output file " << outputFile->GetName() << " successfully opened." << std::endl;
  }
  outputFile->cd();


  ///initialise PU density object

  HGCSSPUenergy puDensity("data/EnergyDensity.dat");


    //////////////////////////////////////////////////
    //////////////////////////////////////////////////
    ///////// positionFit /////////////////////////////
    //////////////////////////////////////////////////
    //////////////////////////////////////////////////
  
  const unsigned nEvts = ((pNevts > lSimTree->GetEntries() || pNevts==0) ? static_cast<unsigned>(lSimTree->GetEntries()) : pNevts) ;
  

  PositionFit lChi2Fit(nSR,residualMax,nLayers,nSiLayers,applyPuMixFix,debug);
  lChi2Fit.initialise(outputFile,"PositionFit",outFolder,geomConv,puDensity);

  //try getting z position from input file, if doesn't exit,
  //perform first loop over simhits to find z positions of layers
  if ((redoStep<2 && !lChi2Fit.getZpositions(versionNumber)) || redoStep>1)
    lChi2Fit.getZpositions(versionNumber,lSimTree,nEvts);
  
  //perform second loop over events to find positions to fit and get energies
  SignalRegion SignalEnergy(outFolder, nLayers, nEvts, geomConv, puDensity,applyPuMixFix,versionNumber);
  SignalEnergy.initialise(outputFile,"Energies");

  //initialise
  bool dofit = redoStep>0 || !SignalEnergy.initialiseFitPositions();
  if (!dofit && redoStep==0) std::cout << " -- Info: fit positions taken from file on disk." << std::endl;
  else std::cout << " -- Info: redoing least square fit." << std::endl;

  if (redoStep>0 || (redoStep==0 && dofit)){
    lChi2Fit.getInitialPositions(lSimTree,lRecTree,nEvts);
    lChi2Fit.finaliseErrorMatrix();
    lChi2Fit.initialiseLeastSquareFit();
  }
  
  //loop on events
  HGCSSEvent * event = 0;
  std::vector<HGCSSSamplingSection> * ssvec = 0;
  std::vector<HGCSSSimHit> * simhitvec = 0;
  std::vector<HGCSSRecoHit> * rechitvec = 0;
  std::vector<HGCSSGenParticle> * genvec = 0;
  unsigned nPuVtx = 0;

  lSimTree->SetBranchAddress("HGCSSEvent",&event);
  lSimTree->SetBranchAddress("HGCSSSamplingSectionVec",&ssvec);
  lSimTree->SetBranchAddress("HGCSSSimHitVec",&simhitvec);
  lSimTree->SetBranchAddress("HGCSSGenParticleVec",&genvec);

  lRecTree->SetBranchAddress("HGCSSRecoHitVec",&rechitvec);
  if (lRecTree->GetBranch("nPuVtx")) lRecTree->SetBranchAddress("nPuVtx",&nPuVtx);
  const unsigned nRemove = 12;
  std::vector<unsigned> lToRemove;
  unsigned list[nRemove] = {25,27,15,1,10,3,18,5,12,7,23,20};

  for (unsigned ievt(0); ievt<nEvts; ++ievt){//loop on entries
    if (debug) std::cout << "... Processing entry: " << ievt << std::endl;
    else if (ievt%50 == 0) std::cout << "... Processing entry: " << ievt << std::endl;

    lSimTree->GetEntry(ievt);
    lRecTree->GetEntry(ievt);
    if (dofit) {
      bool found = lChi2Fit.setTruthInfo(genvec,1);
      if (!found) continue;
      //mask layers in turn
      lToRemove.clear();
      for (unsigned r(0); r<nRemove+1;++r){
	FitResult fit;
	if ( lChi2Fit.performLeastSquareFit(ievt,fit,lToRemove)==0 ){
	  //SignalEnergy.fillEnergies(ievt,(*ssvec),(*simhitvec),(*rechitvec),nPuVtx,fit);
	}
	else std::cout << " -- remove " << r << " Fit failed." << std::endl;
	if (r<nRemove) lToRemove.push_back(list[r]);
      }

    }
    else SignalEnergy.fillEnergies(ievt,(*ssvec),(*simhitvec),(*rechitvec),nPuVtx);

  }//loop on entries

  //finalise

  if (dofit) lChi2Fit.finaliseFit();
  SignalEnergy.finalise();

  outputFile->Write();
  //outputFile->Close();
  
  std::cout << " - End of egammaResolution program." << std::endl;

  return 0;
  

}//main

void PDFSystCalculator_Powheg12Bin()
{
  typedef std::pair<double,double> pairOfDouble;
  set<pairOfDouble> uniqueEventsTrue;

  typedef std::pair<int,int> pairOfInt;
  set<pairOfInt> uniqueEventsReco;

  double mLow = 60.; // Z mass
  double mHigh = 120.; // Z mass
  double ptLow = 0.;
  double ptHigh = 4000.;

  //variable related PDF systematic uncertainty
  //double weightedSelectedEvents[18][53];
  //double weighted2SelectedEvents[18][53];
  double weightedSelectedEvents[12][53];
  double weighted2SelectedEvents[12][53];
  
  double events_central[nptBins] = {0.0};
  double events2_central[nptBins] = {0.0};
  double wa[nptBins] = {0.0};
  double wb[nptBins] = {0.0};
  double wplus[nptBins] = {0.0};
  double wminus[nptBins] = {0.0};

  char tmpName[30];


 // std::ofstream Fout;
 // TString FoutName = "Zpt_12BinWeightedSelectedEvents.txt";
 // Fout.open(FoutName);
 // Fout << fixed << setprecision(3);



  // read input root file
  TChain* tree = new TChain("tree");
  for(int i(1);i<95;i++)
  {
    //sprintf(tmpName,"DYJetsToLL_%d.root",i);
    //sprintf(tmpName,"/d1/scratch/sangilpark/Zpt/CMSSW_5_3_14_patch1/src/TerraNova/NtupleMaker/test/ZpT_Powheg_LowPU/DYJetsToLL_%d.root",i);
    sprintf(tmpName,"/d2/scratch/Storage_Area/ZpT8TeV_Powheg_LowPU/DYJetsToLL_%d.root",i);
    tree->AddFile(tmpName);  // Powheg Ntuple
  }

  TH1D* hPtMCTot = new TH1D("hPtMCTot", "", nptBins, xbins_pt); // no cut
  TH1D* hPtMCAcc = new TH1D("hPtMCAcc", "", nptBins, xbins_pt); // acc.cuts
  TH1D* hPtMCEff = new TH1D("hPtMCEff", "", nptBins, xbins_pt); // acc + all selection cuts


  //Hammid MC truth information pre FSR
  _TrackInfo true1PreFSR, true2PreFSR;
  tree->SetBranchAddress("true1PreFSR", &true1PreFSR);
  tree->SetBranchAddress("true2PreFSR", &true2PreFSR);

  float trueMassPreFSR, truePtPreFSR;
  tree->SetBranchAddress("trueMassPreFSR",	&trueMassPreFSR);
  tree->SetBranchAddress("truePtPreFSR",	&truePtPreFSR);

  vector<float> *weights_CT10;
  TBranch *b_weights_CT10;
  weights_CT10 = 0;
  tree->SetBranchAddress("weights_CT10", &weights_CT10, &b_weights_CT10);


  // Initialize variable
  for(int iBin(0); iBin<nptBins; iBin++)
  {
    {
      for(int j=0; j<53; j++)
      {
	weightedSelectedEvents[iBin][j] = 0;
	weighted2SelectedEvents[iBin][j] = 0;
//	cout << "initialized weightedSelectedEvents : " << weightedSelectedEvents[iBin][j] << endl;
      }
    }
  }

  // Print # of events
  cout << "Loop over the " << tree->GetEntries() << " entries ...\n";

  // Fill histogram
  for(int iEvt(0); iEvt<tree->GetEntries(); iEvt++)
  //for(int iEvt(0); iEvt<20; iEvt++)
  {
    if ( (iEvt % 100000)==0 ) cout << "event " << iEvt << endl;
    tree -> GetEntry(iEvt);

    //nocut start
    pairOfDouble trueMassPtPreFSR(trueMassPreFSR,truePtPreFSR);
    if( !uniqueEventsTrue.insert( trueMassPtPreFSR ).second ) continue;

    if (trueMassPreFSR <  mLow) continue;
    if (trueMassPreFSR > mHigh) continue;

    hPtMCTot -> Fill( truePtPreFSR );

    // acc cut start
    if (true1PreFSR.charge == -999) continue;
//.........这里部分代码省略.........

//---------------------------------------------------------------------------//
//Main Method
//---------------------------------------------------------------------------//
int main(int argc, char** argv)
{	
	TApplication *App = new TApplication("Application",(Int_t*)&argc, argv);
	TCanvas *Canvas = new TCanvas("canvas", "Canvas", 640, 640);

	// Set up ROOT as we require.
	SetupROOT();

	// Get list of files to run over. 
	TString fileName("/storage/epp2/phseaj/exercise/basket_2010b.list");
	std::ifstream inputFile(fileName.Data(), ios::in);

	// Declare a TChain for the TGlobalPID module
	TChain *gRecon = new TChain("ReconDir/Global");
	TChain *gGenVtx = new TChain("TruthDir/Vertices");
	// Check if the file exists.
	if (!inputFile.is_open()){
	std::cout << "ERROR: File prod4 files not found!" << std::endl;
		std::cout << " - This file should contain a list of all data files to be processed." << std::endl;
		return 0;
	}
	else{
		std::string curFileName;

		// Add the input files to the TChains.
		//only doing 10 of the basket files, revert to while to do whole run
	//	while(getline(inputFile,curFileName)){
		for(int l = 0; l<10; l++){
			if(getline(inputFile,curFileName)){
				gRecon->Add(curFileName.c_str());
				gGenVtx->Add(curFileName.c_str());
			}
		}
	}

	std::cout << "Got input file(s)." << std::endl;

	//Setup access to the Recon tree
	int NPIDs(0);  // This variable counts the number of particles per event
	int NVtxFGD1(0), NVtxFGD2(0);
        // Declare a TClonesArray to hold objects of type TGlobalPID
 	TClonesArray *globalPIDs = new TClonesArray("ND::TGlobalReconModule::TGlobalPID",50);
	TClonesArray *VtxFGD1 = new TClonesArray("ND::TTruthVerticesModule::TTruthVertex",50);    
	TClonesArray *VtxFGD2 = new TClonesArray("ND::TTruthVerticesModule::TTruthVertex",50);    
    // Associate the right branch in the TTree to the right local variable
	gRecon->SetBranchAddress("NPIDs",&NPIDs);
    gRecon->SetBranchAddress("PIDs",&globalPIDs);
	gGenVtx->SetBranchAddress("VtxFGD1", &VtxFGD1);
	gGenVtx->SetBranchAddress("NVtxFGD1", &NVtxFGD1);
	gGenVtx->SetBranchAddress("VtxFGD2", &VtxFGD2);
	gGenVtx->SetBranchAddress("NVtxFGD2", &NVtxFGD2);
	//check that truthdir and recon have the same number of entries
	if(gRecon->GetEntries() != gGenVtx->GetEntries()) 
		cout<<"not equal entries, probably wrong"<<endl;
	// Loop over the entries in the TChain.

	//========================================================
	//			Declare Graphs n stuff here
	//========================================================

	//adding tclones arrays for use with detectors
	Int_t NTPCs;
	TClonesArray *TPC;

	Int_t NFDGs;
	TClonesArray *FDG;

	Int_t NECALs;
	TClonesArray *ECAL;

	Int_t NPODs;
	TClonesArray *POD;

	Int_t NSMRDs;
	TClonesArray *SMRD;

	//adding a 2d graph general purpose, change titles each time!
	TH1D *graph1 = new TH1D("graph1","Momenta of TPC PIDs from FDGs", 100, -50.0 , 50.0);
	
	Int_t ninteract(0), nfgd(0),nfgdqes(0),fgdqesneu(0);
	//========================================================
	//	end		Declare Graphs n stuff here
	//========================================================

	// Loop over the entries in the TChain. (only 1/1000 of whole entries atm)
	for(unsigned int i = 0; i < gRecon->GetEntries()/10; ++i) {
		if((i+1)%10000 == 0) std::cout << "Processing event: " << (i+1) << std::endl;
		//display status every 10,000 th entry

	// Get an entry for the Recon tree
		gRecon->GetEntry(i);
		gGenVtx->GetEntry(i);
		ND::TGlobalReconModule::TGlobalPID *gTrack = NULL;

		//added new loop for truth vertex
		gGenVtx->GetEntry(i);
		
//.........这里部分代码省略.........

void TreeFiller(string inFiles, string outFile, string histName) {

	
	TChain *origFiles = new TChain("treeVars", "treeVars");
	origFiles->Add(inFiles.c_str());
	
	int nEntries = origFiles->GetEntries();
	cout << "The old chain contains " << nEntries << " entries." << endl;
	
	
	float jet1Eta, jet2Eta, deltaY, genTopPt1, genTopPt2, jet1Mass, jet2Mass, jet1MinMass, jet2MinMass, jet1BDisc, jet2BDisc, jet1SubjetMaxBDisc, jet2SubjetMaxBDisc,
			jet1tau32, jet2tau32, jet1Pt, jet2Pt, jetPtForMistag, mttMass, mttMassPred, mistagWt, mistagWtAll, mistagWtNsubAll, mistagWtNsub, index, cutflow, NNoutput, ptReweight;
	int jet1NSubjets, jet2NSubjets;
		
	
	origFiles->SetBranchAddress("jet1Eta", &jet1Eta);
	origFiles->SetBranchAddress("jet2Eta", &jet2Eta);
	origFiles->SetBranchAddress("deltaY", &deltaY);
	origFiles->SetBranchAddress("jet1Mass", &jet1Mass);
	origFiles->SetBranchAddress("jet2Mass", &jet2Mass);
	origFiles->SetBranchAddress("jet1minMass", &jet1MinMass);
	origFiles->SetBranchAddress("jet2minMass", &jet2MinMass);
	origFiles->SetBranchAddress("jet1Nsubj", &jet1NSubjets);
	origFiles->SetBranchAddress("jet2Nsubj", &jet2NSubjets);
	origFiles->SetBranchAddress("jet1BDisc", &jet1BDisc);
	origFiles->SetBranchAddress("jet2BDisc", &jet2BDisc);
	origFiles->SetBranchAddress("jet1SubjetMaxBDisc", &jet1SubjetMaxBDisc);
	origFiles->SetBranchAddress("jet2SubjetMaxBDisc", &jet2SubjetMaxBDisc);
	origFiles->SetBranchAddress("jet1tau32", &jet1tau32);
	origFiles->SetBranchAddress("jet2tau32", &jet2tau32);
	origFiles->SetBranchAddress("jet1Pt", &jet1Pt);
	origFiles->SetBranchAddress("jet2Pt", &jet2Pt);
	origFiles->SetBranchAddress("jetPtForMistag", &jetPtForMistag);
	origFiles->SetBranchAddress("mttMass", &mttMass);
	origFiles->SetBranchAddress("mttMassPred", &mttMassPred);
	origFiles->SetBranchAddress("index", &index);
	origFiles->SetBranchAddress("cutflow", &cutflow);
	origFiles->SetBranchAddress("NNoutput", &NNoutput);
	origFiles->SetBranchAddress("mistagWt", &mistagWt);
	origFiles->SetBranchAddress("genTopPt1", &genTopPt1);	
	origFiles->SetBranchAddress("genTopPt2", &genTopPt2);	

	
	TFile *newFile = new TFile(outFile.c_str(), "RECREATE");
	TTree *newTree = origFiles->CloneTree(0);
	newTree->Branch("mistagWtNsub", &mistagWtNsub, "mistagWtNsub/F");
	newTree->Branch("mistagWtNsubAll", &mistagWtNsubAll, "mistagWtNsubAll/F");
	newTree->Branch("mistagWtAll", &mistagWtAll, "mistagWtAll/F");
	newTree->Branch("ptReweight", &ptReweight, "ptReweight/F");

	newTree->SetBranchAddress("mistagWtNsub", &mistagWtNsub);
	newTree->SetBranchAddress("misagWtNsubAll", &mistagWtNsubAll);
	newTree->SetBranchAddress("mistagWtAll", &mistagWtAll);
	newTree->SetBranchAddress("ptReweight", &ptReweight);
	
/*	
	TMVA::Reader* reader = new TMVA::Reader();
	reader->AddVariable("jet1Eta", &jet1Eta);
	reader->AddVariable("jet2Eta", &jet2Eta);
	reader->AddVariable("deltaY", &deltaY);
	reader->AddVariable("jet1Mass", &jet1Mass);
	reader->AddVariable("jet2Mass", &jet2Mass);
	reader->AddVariable("jet1BDisc", &jet1BDisc);
	reader->AddVariable("jet2BDisc", &jet2BDisc);
	reader->AddVariable("jet1SubjetMaxBDisc", &jet1SubjetMaxBDisc);
	reader->AddVariable("jet2SubjetMaxBDisc", &jet2SubjetMaxBDisc);
	reader->AddVariable("jet1tau32", &jet1tau32);
	reader->AddVariable("jet2tau32", &jet2tau32);
	reader->AddVariable("jet1Pt", &jet1Pt);
	reader->AddVariable("jet2Pt", &jet2Pt);
	reader->AddVariable("mttMass", &mttMass);

	reader->BookMVA("MLP", "weights/TMVA_tt_Zp_MLP.weightsZp10_cut4.xml");
*/	
	//TFile *mistagFileLow = new TFile("notCSVL_notCSVM_mistag.root");
	//TFile *mistagFileMed = new TFile("CSVL_notCSVM_mistag.root");
	//TFile *mistagFileHi = new TFile("CSVM_mistag.root");
	TFile *mistagFileLow = new TFile("Jan21_mistag.root");//data_LowBscore_mistag_Dec16.root");
	TFile *mistagFileMed = new TFile("Jan21_mistag.root");//data_MedBscore_mistag_Dec16.root");
	TFile *mistagFileHi = new TFile("Jan21_mistag.root");//data_HiBscore_mistag_Dec16.root");
	TFile *mistagFile = new TFile("Jan21_mistag.root");//data_AllBscore_mistag_Dec16.root");
	histName = "MISTAG_RATE_SUB_TTBAR_Inclusive";
	TH1F *mistagRateHistAll = (TH1F *) mistagFile->Get( histName.c_str() )->Clone();	
	TH1F *mistagRateHistLow = (TH1F *) mistagFileLow->Get( histName.c_str() )->Clone();	
	TH1F *mistagRateHistMed = (TH1F *) mistagFileMed->Get( histName.c_str() )->Clone();	
	TH1F *mistagRateHistHi = (TH1F *) mistagFileHi->Get( histName.c_str() )->Clone();	
	string histName2 = "MISTAG_RATE_SUB_TTBAR_InclNsub";
	TH1F *mistagRateHistNSAll = (TH1F *) mistagFile->Get( histName2.c_str() )->Clone();	
	TH1F *mistagRateHistNSLow = (TH1F *) mistagFileLow->Get( histName2.c_str() )->Clone();	
	TH1F *mistagRateHistNSMed = (TH1F *) mistagFileMed->Get( histName2.c_str() )->Clone();	
	TH1F *mistagRateHistNSHi = (TH1F *) mistagFileHi->Get( histName2.c_str() )->Clone();	
	cout << histName << endl;
	cout << "Entries " << mistagRateHistAll->Integral() << endl;	
	cout << "Entries2 " << mistagRateHistNSHi->Integral() << endl;
	
	for (int i = 0; i < origFiles->GetEntries(); i++){
	
		origFiles->GetEntry(i);
		if  (i % 1000000 == 0) cout << 100*(float(i) / float(nEntries)) << " Percent Complete." << endl;
		mistagWt = 0.000;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  CryEnergyDeposited  = new std::vector<float>  [numOfCry];
  pCryEnergyDeposited = new std::vector<float>* [numOfCry];
  PosXEnDep           = new std::vector<float>  [numOfCry];
  pPosXEnDep          = new std::vector<float>* [numOfCry];
  PosYEnDep           = new std::vector<float>  [numOfCry];
  pPosYEnDep          = new std::vector<float>* [numOfCry];
  PosZEnDep           = new std::vector<float>  [numOfCry];
  pPosZEnDep          = new std::vector<float>* [numOfCry];
  
//   short RunDetectorHit[16];
  
  
  std::vector<float> **pEdep;
  std::vector<float> **px;
  std::vector<float> **py;
  std::vector<float> **pz;
  
  pEdep = new std::vector<float>* [numOfCry];
  px    = new std::vector<float>* [numOfCry];
  py    = new std::vector<float>* [numOfCry];
  pz    = new std::vector<float>* [numOfCry];
  
  for (int i = 0 ; i < numOfCry ; i++)
  {
    pEdep[i] = 0; 
    px[i] = 0;
    py[i] = 0;
    pz[i] = 0;
  }

  Short_t  *detector;
  detector = new Short_t [numOfCh];
  
  tree->SetBranchAddress("Seed",&Seed);
  tree->SetBranchAddress("Run",&Run);
  tree->SetBranchAddress("Event",&Event);
  tree->SetBranchAddress("totalEnergyDeposited",&totalEnergyDeposited);
  tree->SetBranchAddress("NumOptPhotons",&NumOptPhotons);
  tree->SetBranchAddress("NumCherenkovPhotons",&NumCherenkovPhotons);
  
  for (int i = 0 ; i < numOfCry ; i++)
  {
    std::stringstream snames;
    snames << "cry" << i;
    tree->SetBranchAddress(snames.str().c_str(),&pEdep[i]);
    snames.str("");
    snames<< "cry" << i << "PosXEnDep";    
    tree->SetBranchAddress(snames.str().c_str(),&px[i]);
    snames.str("");
    snames<< "cry" << i << "PosYEnDep";
    tree->SetBranchAddress(snames.str().c_str(),&py[i]);
    snames.str("");
    snames<< "cry" << i << "PosZEnDep";
    tree->SetBranchAddress(snames.str().c_str(),&pz[i]);


  }
  for (int i = 0 ; i < numOfCh ; i++)
  {
    std::stringstream snames;
    snames << "detector" << i;
    tree->SetBranchAddress(snames.str().c_str(),&detector[i]);
  }
  
  
  

void process()
{
   Int_t raw[512];   // buffer for input signal and bkg trees
   // buffers for output trees
   Int_t sig[512];
   Int_t cmsig[512];
   Int_t cm[16];

   // pedestal
   const char* fbkg_name = "Raw_Data_FZ320P_05_MSSD_2_250V_K237_Pedestal.dat-events.root";
   TFile* fbkg = TFile::Open(fbkg_name);
   if (!fbkg) cout<< "File not found: " << fbkg <<endl<<exitl;

   TTree* tree = (TTree*) fbkg->Get("etree");
   tree->SetBranchAddress("raw", &raw);

   TH2* h2d = (TH2*) fbkg->Get("h2d");
   new TCanvas;
   h2d->Draw();

   TProfile* profile = (TProfile*) fbkg->Get("profile"); 
   //new TCanvas;
   //profile->Draw();

   Int_t pedestal[512];
   for (int i=0; i<512; i++) {
      // pedestal[i] = profile->GetBinContent(i+1) - 0.5;
      // pedestal[i] = profile->GetBinContent(i+1) + 0.5;
      pedestal[i] = profile->GetBinContent(i+1);
   }

   // cout << "\nPedestals for every channel\n" << endl;
   // for (int i=0; i<512; i++) {
   //    cout << pedestal[i] << " ";
   //    if (i>0 && (i+1)%128==0)
   //       cout << endl;
   // }

   cout<< "processing bkg" <<endl;

   // output file with tree
   const char* obfname = "FZ320P_05_MSSD_2-bkg.root";
   TFile* obfile = TFile::Open(obfname, "recreate");

   TTree* btree = new TTree("btree", "btree");
   btree->Branch("sig", &sig, "sig[512]/I");
   btree->Branch("cmsig", &cmsig, "cmsig[512]/I");
   btree->Branch("cm", &cm, "cm[16]/I");
   btree->SetMarkerStyle(6);
   btree->SetMarkerColor(2);

   for (int jentry=0; jentry<tree->GetEntries(); ++jentry)
   {
      tree->GetEvent(jentry);
      // sig
      for (int i=0; i<512; ++i) {
         sig[i] = raw[i] - pedestal[i];
      }
      // calc common mode
      Int_t group32[32];
      Int_t index32[32];
      for (int igroup=0; igroup<16; ++igroup) {
         for (int istrip=0; istrip<32; ++istrip)   // istrip is number inside group of 32
         {
            group32[istrip] = sig[igroup*32 + istrip];
         }
         // sort array group32 in ascending order
         TMath::Sort(32, group32, index32, kFALSE);
         Int_t median = group32[index32[14]];
         cm[igroup] = median;
      }
      // subtract common mode
      for (int istrip=0; istrip<512; ++istrip) {
         Int_t igroup = istrip/32;
         cmsig[istrip] = sig[istrip] - cm[igroup];
      }
      // Fill sig, cmsig, cm
      btree->Fill();
   }
   obfile->Write();

   /////////////////////////////////////////////////
   //
   // signal tree
   //
   /////////////////////////////////////////////////

   cout<< "processing signal" <<endl;
   
   TChain* chain = new TChain("etree");
   chain->Add("Raw_Data_FZ320P_05_MSSD_250V_K237_Position_1.dat-events.root");
   chain->Add("Raw_Data_FZ320P_05_MSSD_250V_K237_Position_2.dat-events.root");
   chain->Add("Raw_Data_FZ320P_05_MSSD_2_250V_K237_Position_3.dat-events.root");
   chain->Add("Raw_Data_FZ320P_05_MSSD_2_250V_K237_Position_4.dat-events.root");
   chain->SetBranchAddress("raw", &raw);

   // output file with tree
   const char* osfname = "FZ320P_05_MSSD_2-signal.root";
   TFile* osfile = TFile::Open(osfname, "recreate");

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

//.........这里部分代码省略.........
  else if (theSample == "WgammaStar") {
    tree->Add(filesPath + "/" + systematic + "/" + "latino_082_WGstarToElNuMad.root");
    tree->Add(filesPath + "/" + systematic + "/" + "latino_083_WGstarToMuNuMad.root");
    tree->Add(filesPath + "/" + systematic + "/" + "latino_084_WGstarToTauNuMad.root");
  }
  else if (theSample == "HWW125") { 
    tree->Add(filesPath + "/" + systematic + "/" + "latino_1125_ggToH125toWWTo2LAndTau2Nu.root");
  }
  else if (theSample =="GamGamWW"){
    tree->Add(filesPath + "/" + systematic + "/" + "latino_008_GamGamWW.root");    
  }
  else if (theSample == "WJets") {
    tree->Add(filesPath + "/" + systematic + "/" + "latino_080_WJetsToLNuMad.root");
    //CHECK PATH: tree->Add("/gpfs/csic_projects/tier3data/LatinosSkims/ReducedTrees/R53X_S1_V08_S2_V09_S3_V13/MC_LooseLoose/4L/latino_080_WJetsToLNuMad.root");
  }
  else if (theSample == "VVV") {
    tree->Add(filesPath + "/" + systematic + "/" + "latino_088_WWGJets.root");
    tree->Add(filesPath + "/" + systematic + "/" + "latino_089_WZZJets.root");
    tree->Add(filesPath + "/" + systematic + "/" + "latino_090_ZZZJets.root");
    tree->Add(filesPath + "/" + systematic + "/" + "latino_091_WWZJets.root");
    tree->Add(filesPath + "/" + systematic + "/" + "latino_092_WWWJets.root");
    tree->Add(filesPath + "/" + systematic + "/" + "latino_093_TTWJets.root");
    tree->Add(filesPath + "/" + systematic + "/" + "latino_094_TTZJets.root");
    tree->Add(filesPath + "/" + systematic + "/" + "latino_095_TTWWJets.root");
    tree->Add(filesPath + "/" + systematic + "/" + "latino_096_TTGJets.root");
  }
  else {
    return;
  }


  // Declaration of leaf types
  //----------------------------------------------------------------------------
  Float_t baseW;        tree->SetBranchAddress("baseW"       , &baseW);
  Float_t channel;      tree->SetBranchAddress("channel"     , &channel);
  Float_t chmet;        tree->SetBranchAddress("chmet"       , &chmet);
  Float_t dataset;      tree->SetBranchAddress("dataset"     , &dataset);
  Float_t dphill;       tree->SetBranchAddress("dphill"      , &dphill);
  Float_t dphilljet;    tree->SetBranchAddress("dphilljet"   , &dphilljet);
  Float_t dphilljetjet; tree->SetBranchAddress("dphilljetjet", &dphilljetjet);
  Float_t drll;         tree->SetBranchAddress("drll"        , &drll);
  Float_t effW;         tree->SetBranchAddress("effW"        , &effW);
  Float_t jeteta1;      tree->SetBranchAddress("jeteta1"     , &jeteta1);
  Float_t jeteta2;      tree->SetBranchAddress("jeteta2"     , &jeteta2);
  Float_t jeteta3;      tree->SetBranchAddress("jeteta3"     , &jeteta3);
  Float_t jetpt1;       tree->SetBranchAddress("jetpt1"      , &jetpt1);
  Float_t jetpt2;       tree->SetBranchAddress("jetpt2"      , &jetpt2);
  Float_t jetpt3;       tree->SetBranchAddress("jetpt3"      , &jetpt3);
  Float_t jettche1;     tree->SetBranchAddress("jettche1"    , &jettche1);
  Float_t jettche2;     tree->SetBranchAddress("jettche2"    , &jettche2);
  Float_t mctruth;      tree->SetBranchAddress("mctruth"     , &mctruth);
  Float_t mll;          tree->SetBranchAddress("mll"         , &mll);
  Float_t mpmet;        tree->SetBranchAddress("mpmet"       , &mpmet); 
  Float_t mth;          tree->SetBranchAddress("mth"         , &mth);
  Float_t nbjet;        tree->SetBranchAddress("nbjet"       , &nbjet);
  Float_t nbjettche;    tree->SetBranchAddress("nbjettche"   , &nbjettche);
  Float_t nextra;       tree->SetBranchAddress("nextra"      , &nextra);
  Float_t njet;         tree->SetBranchAddress("njet"        , &njet);
  Float_t nvtx;         tree->SetBranchAddress("nvtx"        , &nvtx);
  Float_t pchmet;       tree->SetBranchAddress("pchmet"      , &pchmet);
  Float_t pfmet;        tree->SetBranchAddress("pfmet"       , &pfmet);
  Float_t ppfmet;       tree->SetBranchAddress("ppfmet"      , &ppfmet);
  Float_t pt1;          tree->SetBranchAddress("pt1"         , &pt1);
  Float_t pt2;          tree->SetBranchAddress("pt2"         , &pt2);
  Float_t eta1;          tree->SetBranchAddress("eta1"         , &eta1);
  Float_t eta2;          tree->SetBranchAddress("eta2"         , &eta2);

void draw() {

  TChain *chain = new TChain("Events");

  chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_1.root");
  chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_10.root");
  chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_100.root");
  chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_101.root");
  chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_102.root");
  chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_103.root");
  chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_104.root");
  chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_105.root");
  chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_106.root");
  chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_107.root");
  chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_108.root");
  chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_109.root");
  chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_11.root");
  chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_110.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_111.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_112.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_113.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_114.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_115.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_116.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_117.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_118.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_119.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_12.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_120.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_121.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_122.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_123.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_124.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_125.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_126.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_127.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_128.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_129.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_13.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_130.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_131.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_132.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_133.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_134.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_135.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_136.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_137.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_138.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_139.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_14.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_140.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_141.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_142.root");
  //chain->AddFile("/afs/cern.ch/work/j/jlawhorn/HighPtJet80/HCALTree_143.root");
  
  cout << chain->GetEntries() << endl;

  Double_t pulse[10], ped[10];

  chain->SetBranchAddress("pulse", &pulse);
  chain->SetBranchAddress("ped", &ped);

  vector<Double_t> ucut;
  vector<Double_t> lcut;

  const int nbin=11;

  for (int i=0; i<nbin; i++) {
    lcut.push_back(25+i*25);
    ucut.push_back(50+i*25);
  }
  
  Double_t avgPulse[nbin][10], avgPulse2[nbin][10], adcUnc[nbin][10], totUnc[nbin][10]; 
  for (UInt_t i=0; i<nbin; i++) {
    for (UInt_t j=0; j<10; j++) {avgPulse[i][j]=0; avgPulse2[i][j]=0; adcUnc[i][j]=0; totUnc[i][j]=0;}
  }

  Int_t nPulses[nbin]={0};

  TFile *outf = new TFile("test_hist.root", "recreate");
  TH1D *hQ = new TH1D("hq", "hq",50, 0, 400);

  for (UInt_t i=0; i<chain->GetEntries(); i++) {
  //for (UInt_t i=0; i<1000; i++) {
    chain->GetEntry(i);

    Double_t sumQ=0;
    for (UInt_t j=0; j<10; j++) { sumQ+=pulse[j];}

    hQ->Fill(sumQ);

    if (sumQ<25 || sumQ>300) continue;
    Int_t ibin= floor( sumQ/25-1 );
    //cout << sumQ << ", " << ibin << endl;

    nPulses[ibin]++;

    for (UInt_t j=0; j<10; j++) { 
      avgPulse[ibin][j]+=pulse[j]/sumQ;
      avgPulse2[ibin][j]+=pulse[j]*pulse[j]/sumQ;
//.........这里部分代码省略.........

void createPlot(std::vector<TString> samples, std::vector<TString> files, std::vector<TString> legend)
{

  TString y_title = "Number of Entries";
  const int nHist = files.size(); // number of files

  // Declare the histograms to be saved 
  TH1F *h_mll[nHist];
  TH1F *h_dilpt[nHist];
  TH1F *h_dphill[nHist];
  TH1F *h_leadleppt[nHist];
  TH1F *h_trailleppt[nHist];
  TH1F *h_met[nHist];
  TH1F *h_mt[nHist];
  
  // Get the histograms from the ntuples
  for (int i=0;i<nHist;i++) {
    TString treeName = "angles";
    TChain *chain = new TChain(treeName);
    chain->Add(files[i]);
    assert(chain);
    
    // declare histograms  to fill
    Color_t color = kBlack;
    TString sampleName = samples[i];
    if ( sampleName.Contains("SMHiggs",TString::kExact )) color = kBlue;
    if ( sampleName.Contains("PSHiggs", TString::kExact)) color = kMagenta;
    if ( sampleName.Contains("T",  TString::kExact)) color = kRed;
    if ( sampleName.Contains("THiggs",  TString::kExact)) color = kRed;


    
    // define the histograms to plot
    
    // dilmass 
    h_mll[i] = new TH1F(TString("HWW_"+sampleName+"_hdilmass"), TString("HWW_"+sampleName+"_hdilmass"), 20, 0, 200);
    h_mll[i]->SetLineColor(color);
    h_mll[i]->SetMarkerColor(color);
    
    // leading lepton pT
    h_leadleppt[i] = new TH1F(TString("HWW_"+sampleName+"_hleadleppt"), TString("HWW_"+sampleName+"_hleadleppt"), 20, 0, 100);
    h_leadleppt[i]->SetLineColor(color);
    h_leadleppt[i]->SetMarkerColor(color);
   
    // trailing lepton pT
    h_trailleppt[i] = new TH1F(TString("HWW_"+sampleName+"_htrailleppt"), TString("HWW_"+sampleName+"_htrailleppt"), 20, 0, 100);
    h_trailleppt[i]->SetLineColor(color);
    h_trailleppt[i]->SetMarkerColor(color);
    
    // MET
    h_met[i] = new TH1F(TString("HWW_"+sampleName+"_hmet"), TString("HWW_"+sampleName+"_hmet"), 20, 0, 100);
    h_met[i]->SetLineColor(color);
    h_met[i]->SetMarkerColor(color);
    
    
    // dilepton pT
    h_dilpt[i] = new TH1F(TString("HWW_"+sampleName+"_hdilpt"), TString("HWW_"+sampleName+"_hdilpt"), 30, 20, 100);
    h_dilpt[i]->SetLineColor(color);
    h_dilpt[i]->SetMarkerColor(color);
    
    // deltaphi (ll)
    h_dphill[i] = new TH1F(TString("HWW_"+sampleName+"_hdphi"), TString("HWW_"+sampleName+"_hdphi"), 18, 0, 180.0);
    h_dphill[i]->SetLineColor(color);
    h_dphill[i]->SetMarkerColor(color);
    
    // transverse mass
    h_mt[i] = new TH1F(TString("HWW_"+sampleName+"_hmt"), TString("HWW_"+sampleName+"_hmt"), 20, 0, 200);
    h_mt[i]->SetLineColor(color);
    h_mt[i]->SetMarkerColor(color);
    

    std::cout  << "Processing " << chain->GetEntries() << " entries. \n";
    int nEntries =  chain->GetEntries() ;
    int nSelected = 0;
    
    // mcfm variables to be used
    double mll_ = 0.0;
    double leadleppt_ = 0.0;
    double trailleppt_ = 0.0;
    double leadlepeta_ = 0.0;
    double traillepeta_ = 0.0;
    double dphill_ = 0.0;
    double met_ = 0.0;
    double mt_ = 0.0;
    double dilpt_ = 0.0;
    double wt_ = 1.0;

    if (chain->GetBranchStatus("mll"))
      chain->SetBranchAddress("mll", &mll_);

    if (chain->GetBranchStatus("leadleppt"))
      chain->SetBranchAddress("leadleppt", &leadleppt_);

    if (chain->GetBranchStatus("trailleppt"))
      chain->SetBranchAddress("trailleppt", &trailleppt_);

    if (chain->GetBranchStatus("dphill"))
      chain->SetBranchAddress("dphill", &dphill_);

    if (chain->GetBranchStatus("met"))
//.........这里部分代码省略.........

void _tag()
{
	TCanvas * c1 = new TCanvas("c1", "The 3d view",0,0,1000,500);
	c1->Divide(2,1);
	TChain* T2 = new TChain("Stats");
	T2->Add("TrashMCTest.root");

	TChain* T = new TChain("Primaries");
	T->Add("VertexRestorer.root");
	
	float _primeCostheta[MAXN];
	int _tag = 0;
	int _primariesTotal = 0;

	T2->SetBranchAddress("tag", &_tag);
	T->SetBranchAddress("primariesTotal", &_primariesTotal);
	T->SetBranchAddress("primeCostheta", &_primeCostheta);

	int mTotalNumberOfEvents = T2->GetEntries();
	int nbins = 30;
	int maxd = 1;
	TH1F * probhist = new TH1F("Pw","# of tracks comparison",nbins,-1,maxd);
	TH1F * taghist = new TH1F("Pw","# of tracks comparison",nbins,-1,maxd);
	for (int i = 0; i < mTotalNumberOfEvents; i++) 
	{
		T2->GetEntry(i);
		T->GetEntry(i);
		if (_tag == 0) 
		{
			for (int j = 0; j < _primariesTotal; j++) 
			{
				probhist->Fill(_primeCostheta[j]);
			}
		}
		else 
		{
			for (int j = 0; j < _primariesTotal; j++) 
			{
				taghist->Fill(_primeCostheta[j]);
			}
		}
	
	}
	c1->cd(1);
	THStack * stack = new THStack("d",";cos");
	stack->Add(probhist);
	stack->Add(taghist);
	probhist->SetLineWidth(3);
	//probhist->SetFillColor(kBlue);
	taghist->SetLineWidth(3);
	taghist->SetLineColor(kYellow+1);
	//taghist->SetFillColor(kYellow);
	probhist->Draw();
	taghist->Draw("same");
	stack->Draw();
	c1->cd(2);
	TFile * file = TFile::Open("TrashRecoTest.root");
	file->cd();
	TTree * tree = Jets;
	TH1F * bhist = new TH1F("bhist","# of tracks comparison",nbins,-1,maxd);
	TH1F * bbarhist = new TH1F("bbarhist","# of tracks comparison",nbins,-1,maxd);
	TH1F * whist = new TH1F("whist",";cos#theta",nbins,-1,maxd);
	tree->Project("bhist", "costhetaJetParticles", "costhetaJetParticles > -1 && mcpdg < 0");
	tree->Project("bbarhist", "costhetaJetParticles", "costhetaJetParticles > -1 && mcpdg > 0");
	tree->Project("whist", "costhetaJetParticles", "costhetaJetParticles > -1 && mcpdg == 0");
	whist->SetLineWidth(3);
	whist->SetLineColor(kBlue);
	whist->SetMinimum(0);
	whist->Draw("same");
	bhist->SetLineWidth(3);
	bhist->SetLineColor(kGreen);
	bhist->Draw("same");
	bbarhist->SetLineWidth(3);
	bbarhist->SetLineColor(kRed);
	bbarhist->Draw("same");
	TLegend *legendMean = new TLegend(0.17,0.7,0.5,0.92,NULL,"brNDC");
	legendMean->SetFillColor(kWhite);
	legendMean->SetBorderSize(0);
	legendMean->AddEntry(bhist,"b-jet","fp");
	legendMean->AddEntry(bbarhist,"#bar{b}-jet","fp");
	legendMean->AddEntry(whist,"w-jets","fp");
	legendMean->Draw();
	gPad->Modified();

}

void fitSingleMass( const std::string& basedir, float mass, const std::string&  width, TGraphErrors* gr_mean, TGraphErrors* gr_sigma, TGraphErrors* gr_width, TGraphErrors* gr_alpha1, TGraphErrors* gr_n1, TGraphErrors* gr_alpha2, TGraphErrors* gr_n2 ) {


  std::string outdir = "genSignalShapes";
  system( Form("mkdir -p %s", outdir.c_str()) );


  std::string dataset( Form( "GluGluSpin0ToZGamma_ZToLL_W_%s_M_%.0f_TuneCUEP8M1_13TeV_pythia8", width.c_str(), mass ) );
  std::cout << "-> Starting: " << dataset << std::endl;

  system( Form("ls %s/%s/crab_%s/*/0000/genAna_1.root >> toBeAdded.txt", basedir.c_str(), dataset.c_str(), dataset.c_str() ) );
  TChain* tree = new TChain("mt2");
  ifstream ifs("toBeAdded.txt");
  while( ifs.good() ) {
    std::string fileName;
    ifs >> fileName;
    TString fileName_tstr(fileName);
    if( !fileName_tstr.Contains("pnfs") ) continue;
    tree->Add(Form("$DCAP/%s/mt2", fileName.c_str()) );
  }
  system( "rm toBeAdded.txt" );

  if( tree->GetEntries()==0 ) return;

  int ngenPart;
  tree->SetBranchAddress( "ngenPart", &ngenPart );
  float genPart_pt[100];
  tree->SetBranchAddress( "genPart_pt", genPart_pt );
  float genPart_eta[100];
  tree->SetBranchAddress( "genPart_eta", genPart_eta );
  float genPart_phi[100];
  tree->SetBranchAddress( "genPart_phi", genPart_phi );
  float genPart_mass[100];
  tree->SetBranchAddress( "genPart_mass", genPart_mass );
  int genPart_pdgId[100];
  tree->SetBranchAddress( "genPart_pdgId", genPart_pdgId );
  int genPart_motherId[100];
  tree->SetBranchAddress( "genPart_motherId", genPart_motherId );
  int genPart_status[100];
  tree->SetBranchAddress( "genPart_status", genPart_status );



  RooRealVar* x = new RooRealVar("boss_mass", "boss_mass", mass, 0.5*mass, 1.5*mass );

  RooDataSet* data = new RooDataSet( "data", "data", RooArgSet(*x) );



  int nentries = tree->GetEntries();


  for( int iEntry = 0; iEntry<nentries; ++iEntry ) {

    if( iEntry % 25000 == 0 ) std::cout << "  Entry: " << iEntry << " / " << nentries << std::endl;

    tree->GetEntry(iEntry);

    TLorentzVector leptPlus;
    TLorentzVector leptMinus;
    TLorentzVector photon;
    bool foundLeptPlus = false;
    bool foundLeptMinus = false;
    bool foundPhoton = false;
    bool tauEvent = false;

    for( int iPart=0; iPart<ngenPart; ++iPart ) {

      if( genPart_status[iPart]!=1 ) continue;

      if( abs(genPart_pdgId[iPart])==15 ) {
        tauEvent = true;
        break;
      }

      if( (genPart_pdgId[iPart]==+11 || genPart_pdgId[iPart]==+13) && genPart_motherId[iPart]==23 ) {
        leptMinus.SetPtEtaPhiM( genPart_pt[iPart], genPart_eta[iPart], genPart_phi[iPart], genPart_mass[iPart] );
        foundLeptMinus = true;
      }
      if( (genPart_pdgId[iPart]==-11 || genPart_pdgId[iPart]==-13) && genPart_motherId[iPart]==23 ) {
        leptPlus.SetPtEtaPhiM( genPart_pt[iPart], genPart_eta[iPart], genPart_phi[iPart], genPart_mass[iPart] );
        foundLeptPlus = true;
      }
      if( genPart_pdgId[iPart]==22 && genPart_motherId[iPart]==25 ) {
        photon.SetPtEtaPhiM( genPart_pt[iPart], genPart_eta[iPart], genPart_phi[iPart], genPart_mass[iPart] );
        foundPhoton = true;
      }

    } // for genparts

    if( tauEvent ) continue;
    if( !foundLeptPlus || !foundLeptMinus || !foundPhoton ) continue;


    if( photon.Pt() < 40. ) continue;
    float ptMax = TMath::Max( leptPlus.Pt(), leptMinus.Pt() );
    float ptMin = TMath::Min( leptPlus.Pt(), leptMinus.Pt() );
    if( ptMax<25. ) continue;
    if( ptMin<20. ) continue;
    if( fabs( photon.Eta() ) > 2.5 ) continue;
//.........这里部分代码省略.........

void fit_JVBF(int erg_tev=8, float mSample=0){
  float mPOLE=mSample;
  if (mPOLE<=0) mPOLE=125.6;
  float wPOLE=4.15e-3;
  char TREE_NAME[] = "TestTree";

  TString cinput;
  if (mSample>0) cinput = Form("HZZ4lTree_jvbfMELA_H%.0f_%iTeV.root", mSample, erg_tev);
  else cinput = Form("HZZ4lTree_jvbfMELA_HAll_%iTeV.root", erg_tev);

  TFile* foutput;
  if (mSample>0) foutput = new TFile(Form("jvbfMELA_Fits_H%.0f_%iTeV.root", mSample, erg_tev), "recreate");
  else foutput = new TFile(Form("jvbfMELA_fits_wide_%iTeV.root", erg_tev), "recreate");

  TLorentzVector nullFourVector(0, 0, 0, 0);

  float MC_weight_noxsec;

  float pjvbf_VAJHU;
  float pjvbf_VAJHU_first;
  float pjvbf_VAJHU_second;

  float jet1Pt, jet2Pt;
  float jet1Eta, jet2Eta;
  float jet1Phi, jet2Phi;
  float jet1E, jet2E;
  float jet1Pt_Fake, jet2Pt_Fake;
  float jet1Eta_Fake, jet2Eta_Fake;
  float jet1Phi_Fake, jet2Phi_Fake;
  float jet1E_Fake, jet2E_Fake;
  float jet1px, jet1py, jet1pz;
  float jet2px, jet2py, jet2pz;
  float ZZPx, ZZPy, ZZPz, ZZE, dR;
  short NJets30;

  float ZZMass;

  TChain* tree = new TChain(TREE_NAME);
  tree->Add(cinput);
  tree->SetBranchAddress("MC_weight_noxsec", &MC_weight_noxsec);
  tree->SetBranchAddress("ZZMass", &ZZMass);
  tree->SetBranchAddress("pjvbf_VAJHU", &pjvbf_VAJHU);
  tree->SetBranchAddress("pjvbf_VAJHU_first", &pjvbf_VAJHU_first);
  tree->SetBranchAddress("pjvbf_VAJHU_second", &pjvbf_VAJHU_second);
  tree->SetBranchAddress("NJets30", &NJets30);
  tree->SetBranchAddress("jet1Pt", &jet1Pt);
  tree->SetBranchAddress("jet1Eta", &jet1Eta);
  tree->SetBranchAddress("jet1Phi", &jet1Phi);
  tree->SetBranchAddress("jet1E", &jet1E);
  tree->SetBranchAddress("jet2Pt", &jet2Pt);
  tree->SetBranchAddress("jet2Eta", &jet2Eta);
  tree->SetBranchAddress("jet2Phi", &jet2Phi);
  tree->SetBranchAddress("jet2E", &jet2E);
  tree->SetBranchAddress("jet1_Fake_Pt", &jet1Pt_Fake);
  tree->SetBranchAddress("jet1_Fake_Eta", &jet1Eta_Fake);
  tree->SetBranchAddress("jet1_Fake_Phi", &jet1Phi_Fake);
  tree->SetBranchAddress("jet1_Fake_E", &jet1E_Fake);
  tree->SetBranchAddress("jet2_Fake_Pt", &jet2Pt_Fake);
  tree->SetBranchAddress("jet2_Fake_Eta", &jet2Eta_Fake);
  tree->SetBranchAddress("jet2_Fake_Phi", &jet2Phi_Fake);
  tree->SetBranchAddress("jet2_Fake_E", &jet2E_Fake);


  const int nbinsx=6;
  double bins_mzz[nbinsx+1]={70, 140, 190, 240, 340, 520, 20000};
  double sum_w[nbinsx]={ 0 };
  double sum_mzz[nbinsx][2]={ { 0 } };
  const int nbins_eta = 8;
  double bins_eta[nbins_eta+1]={ 0, 0.75, 1.5, 2.25, 3, 4, 5.5, 10.5, 12 };

  TProfile* hProb[nbinsx];
  TH1D* hProb_all;
  TProfile* hEta[nbinsx+1];
  TH1D* hDist[nbinsx+1];
  for (int bin=0; bin<nbinsx+1; bin++){
    if (bin<nbinsx) hProb[bin] = new TProfile(Form("njets1_pjvbf_bin_%i", bin+1), Form("m_{4l}: [%.0f, %.0f] GeV", bins_mzz[bin], bins_mzz[bin+1]), nbins_eta, bins_eta);
    else hProb_all = new TH1D("njets1_pjvbf_all", Form("m_{4l}: [%.0f, %.0f] GeV", bins_mzz[0], bins_mzz[nbinsx]), nbins_eta, bins_eta);
    if (bin<nbinsx) hEta[bin] = new TProfile(Form("njets1_eta_bin_%i", bin+1), Form("m_{4l}: [%.0f, %.0f] GeV", bins_mzz[bin], bins_mzz[bin+1]), nbins_eta, bins_eta);
    else hEta[bin] = new TProfile("njets1_eta_all", Form("m_{4l}: [%.0f, %.0f] GeV", bins_mzz[0], bins_mzz[nbinsx]), nbins_eta, bins_eta);
    if (bin<nbinsx) hDist[bin] = new TH1D(Form("njets1_dist_bin_%i", bin+1), Form("m_{4l}: [%.0f, %.0f] GeV", bins_mzz[bin], bins_mzz[bin+1]), nbins_eta, bins_eta);
    else hDist[bin] = new TH1D("njets1_dist_all", Form("m_{4l}: [%.0f, %.0f] GeV", bins_mzz[0], bins_mzz[nbinsx]), nbins_eta, bins_eta);
    if (bin<nbinsx) hProb[bin]->Sumw2();
    else hProb_all->Sumw2();
    hEta[bin]->Sumw2();
    hDist[bin]->Sumw2();
  }

  int nEntries = tree->GetEntries();
  cout << nEntries << endl;
  for (int ev = 0; ev < nEntries; ev++){
    tree->GetEntry(ev);

    int massbin=-1;
    for (int bin=0; bin<nbinsx; bin++){
      if (ZZMass>=bins_mzz[bin] && ZZMass<bins_mzz[bin+1]){
        massbin=bin;
      }
    }
    if (NJets30==1){
      hProb[massbin]->Fill(fabs(jet2Eta_Fake), pjvbf_VAJHU_second, MC_weight_noxsec);
//.........这里部分代码省略.........

void test_JVBF(int erg_tev=8, float mSample=0, bool isggH=false){
  float mPOLE=mSample;
  if (mPOLE<=0) mPOLE=125.6;
  float wPOLE=4.15e-3;
  char TREE_NAME[] = "SelectedTree";

  //	TVar::VerbosityLevel verbosity = TVar::INFO;

  Mela mela(erg_tev, mPOLE);

  TFile* foutput;
  if (!isggH){
    if (mSample>0) foutput = new TFile(Form("HZZ4lTree_jvbfMELA_H%.0f_%iTeV.root", mSample, erg_tev), "recreate");
    else foutput = new TFile(Form("HZZ4lTree_jvbfMELA_HAll_%iTeV.root", erg_tev), "recreate");
  }
  else{
    if (mSample>0) foutput = new TFile(Form("HZZ4lTree_jvbfMELA_ggH%.0f_%iTeV.root", mSample, erg_tev), "recreate");
    else foutput = new TFile(Form("HZZ4lTree_jvbfMELA_ggHAll_%iTeV.root", erg_tev), "recreate");
  }
  TLorentzVector nullFourVector(0, 0, 0, 0);

  float MC_weight_noxsec;

  float pjvbf_VAJHU;
  float pjvbf_VAJHU_first;
  float pjvbf_VAJHU_second;
  float phj_VAJHU_first;
  float phj_VAJHU_second;
  float pAux_vbf;
  float pAux_vbf_first;
  float pAux_vbf_second;

  float jet1Pt, jet2Pt;
  float jet1Eta, jet2Eta;
  float jet1Phi, jet2Phi;
  float jet1E, jet2E;
  float jet1Pt_Fake, jet2Pt_Fake;
  float jet1Eta_Fake, jet2Eta_Fake;
  float jet1Phi_Fake, jet2Phi_Fake;
  float jet1E_Fake, jet2E_Fake;
  float jet1px, jet1py, jet1pz;
  float jet2px, jet2py, jet2pz;
  float ZZPx, ZZPy, ZZPz, ZZE, dR;
  short NJets30;
  std::vector<double> * JetPt=0;
  std::vector<double> * JetEta=0;
  std::vector<double> * JetPhi=0;
  std::vector<double> * JetMass=0;
  std::vector<double> myJetPt;
  std::vector<double> myJetCosTheta;
  std::vector<double> myJetEta;
  std::vector<double> myJetPhi;
  std::vector<double> myJetMass;
  TBranch* bJetPt=0;
  TBranch* bJetEta=0;
  TBranch* bJetPhi=0;
  TBranch* bJetMass=0;

  float ZZMass, ZZPt, ZZPhi, ZZEta;

  int GenLep1Id, GenLep2Id, GenLep3Id, GenLep4Id;

  TChain* tree = new TChain(TREE_NAME);
  char* user_folder[3]={ "4mu", "4e", "2mu2e" };
  TString cinput_main = "/scratch0/hep/ianderso/CJLST/140519/PRODFSR";
  if (erg_tev==8) cinput_main.Append("_8TeV");

//  TString cinput_main = "/afs/cern.ch/work/u/usarica/HZZ4l-125p6-FullAnalysis/LHC_";
//  cinput_main.Append(Form("%iTeV", erg_tev));

  for (int ff=0; ff<3; ff++){
    if (!isggH){
      if (mSample>0) tree->Add(Form("%s/%s/HZZ4lTree_VBFH%.0f.root", cinput_main.Data(), user_folder[ff], mSample));
      else tree->Add(Form("%s/%s/HZZ4lTree_VBFH*.root", cinput_main.Data(), user_folder[ff]));
    }
    else{
      if (mSample>0) tree->Add(Form("%s/%s/HZZ4lTree_minloH%.0f.root", cinput_main.Data(), user_folder[ff], mSample));
      else tree->Add(Form("%s/%s/HZZ4lTree_minloH*.root", cinput_main.Data(), user_folder[ff]));
    }
  }
  tree->SetBranchAddress("MC_weight_noxsec", &MC_weight_noxsec);
  tree->SetBranchAddress("NJets30", &NJets30);
  tree->SetBranchAddress("JetPt", &JetPt, &bJetPt);
  tree->SetBranchAddress("JetEta", &JetEta, &bJetEta);
  tree->SetBranchAddress("JetPhi", &JetPhi, &bJetPhi);
  tree->SetBranchAddress("JetMass", &JetMass, &bJetMass);
  tree->SetBranchAddress("ZZMass", &ZZMass);
  tree->SetBranchAddress("ZZPt", &ZZPt);
  tree->SetBranchAddress("ZZEta", &ZZEta);
  tree->SetBranchAddress("ZZPhi", &ZZPhi);

  TTree* newtree = new TTree("TestTree", "");
  newtree->Branch("MC_weight_noxsec", &MC_weight_noxsec);
  newtree->Branch("ZZMass", &ZZMass);
  newtree->Branch("pAux_vbf", &pAux_vbf);
  newtree->Branch("pAux_vbf_first", &pAux_vbf_first);
  newtree->Branch("pAux_vbf_second", &pAux_vbf_second);
  newtree->Branch("pjvbf_VAJHU", &pjvbf_VAJHU);
  newtree->Branch("pjvbf_VAJHU_first", &pjvbf_VAJHU_first);
  newtree->Branch("pjvbf_VAJHU_second", &pjvbf_VAJHU_second);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  }
  
  
  std::cout << " -- 2-D histograms: " << std::endl
	    << " -- X: " << nX << " " << minX << " " << maxX << std::endl
	    << " -- Y: " << nY << " " << minY << " " << maxY << std::endl
	    << " -- Z: " << nZ << " " << minZ << " " << maxZ << std::endl
    ;
  outputFile->cd();

  TH1F *p_nRecHits = new TH1F("p_nRecHits","n(RecHits)",
			      1000,0,500000);
  p_nRecHits->StatOverflows();

  TProfile *p_EvsLayer = new TProfile("p_EvsLayer","Average E vs Layer",30,1,30);

  TProfile *p_AveE[nLayers];
  std::ostringstream p_AveE_Name;
  for(unsigned iL(0); iL<nLayers; iL++){
     p_AveE_Name.str("");
     p_AveE_Name << "p_AveE_layer" <<iL+1;
     p_AveE[iL] = new TProfile(p_AveE_Name.str().c_str(),"Average Energy vs eta",40,1.5,3.5);
  }

  TH2F *p_Occupancy[nLayers];
  for(unsigned iL(0); iL<nLayers; iL++){
     p_AveE_Name.str("");
     p_AveE_Name << "p_Occupancy_Layer" <<iL+1;
     p_Occupancy[iL] = new TH2F(p_AveE_Name.str().c_str(),"RecoHits Occupancy",340,-170,170,340,-170,170);
 }

  std::vector<HGCSSRecoHit> * rechitvec = 0;
  
  lRecTree->SetBranchAddress("HGCSSRecoHitVec",&rechitvec);

  const unsigned nEvts = ((pNevts > lRecTree->GetEntries() || pNevts==0) ? static_cast<unsigned>(lRecTree->GetEntries()) : pNevts) ;
  
  std::cout << "- Processing = " << nEvts  << " events out of " << lRecTree->GetEntries() << std::endl;
  
  //Initialise histos
  //necessary to have overflows ?
  gStyle->SetOptStat(0);

  double Z_layer[nLayers];
  for(unsigned iL(0);iL<nLayers;iL++){
    Z_layer[iL]=0;
  }
  
  double CellNumber[nLayers][40];
  for(unsigned iL(0);iL<nLayers;iL++){
     for(unsigned iEta(0);iEta<40;iEta++){
        CellNumber[iL][iEta] =0;
      }
  }

  std::map<unsigned,bool> channelAlive;

  for (unsigned ievt(0); ievt<nEvts; ++ievt){//loop on entries
    if (debug) std::cout << "... Processing entry: " << ievt << std::endl;
    else if (ievt%50 == 0) std::cout << "... Processing entry: " << ievt << std::endl;
    
    lRecTree->GetEntry(ievt);

    double EvsLayer[nLayers][40];
    double EtotRec[nLayers];
    for(unsigned iL(0);iL<nLayers;iL++){