C++ TTree::Branch方法代码示例

int main(int argc, char* argv[]) {

  if (argc < 2) {
    std::cout << "Insufficient arguments" << std::endl;
    std::cout << "Usage : simulate <params file>" << std::endl;
    exit(1);
  }

  // arguments
  std::string paramfile = std::string(argv[1]);


  Experiment* e;


  // read parameters from file and create experiment object
  e = new Experiment(paramfile);
  std::cout << "Simulation parameters" << std::endl;
  std::cout << "From file : " << paramfile << std::endl;
  std::cout << *e << std::endl;

  // setup tree to write out
  TFile f("outfile.root", "RECREATE");
  TTree *tree = new TTree("T", "Experiments from the Stopped Particle MC");
  tree->Branch("Experiment", "Experiment", &e);

  // fills file
  std::ofstream fillsfile;
  fillsfile.open("fills.txt");

  // setup simulation
  Simulator sim;

  // set parameters
  sim.setParameters(e);

  // set observed events
  sim.setupObservedEvents();

  // set up luminosity data
  sim.setupLumi(false, true, true, e->lumiFile);

  // set up lifetime mask
  sim.setupLifetimeMask();

  //  sim.getLhcFills().print(std::cout);


  // simulate everything at once
  unsigned firstFill = e->fills.at(0);
  unsigned lastFill  = e->fills.at(e->fills.size()-1);

  // expected background
  sim.calculateExpectedBG(firstFill, lastFill);
  
  // get observed events
  sim.calculateObservedEvents(firstFill, lastFill);
  
  // run MC
  sim.simulateSignal(firstFill, lastFill);
  
  // calculate limit
  sim.calculateLimit();

  std::cout << std::endl;

  // write results summary
  std::ofstream summaryfile;
  summaryfile.open("summary.txt", std::ios::app);
  e->summary(summaryfile);
  summaryfile.close();

  std::cout << *e << std::endl;

  tree->Fill();
  f.cd();
  tree->Write();

  fillsfile.close();

  return 0;
}

// provide either a root file or text file with root file paths
int createTrainingTree(TString infile, TString outfile)
{
    //static
    // initialize variables for filling struct 
    Float_t fEnergyGeV, fPrimaryZenithDeg, fPrimaryAzimuthDeg, fLength, fWidth, fSize, fDist, fTimeGradient, fLoss;
    triggeredSimEvent trigSimEvent = triggeredSimEvent(); // for writing 

    // for branch, if type is other than float, need to specify types with varname/F, for example
    const char* leafList = "fEnergyGeV:fLength:fWidth:fSize";
    //:fPrimaryZenithDeg:fPrimaryAzimuthDeg:fTimeGradient:fLoss:fDist
    
    // quality cuts 
    int NTubesMin = 5; // default in VEGAS 
    float distanceUpper = 1.38; // 1.43 [degrees]

    
    // create output file for writing new branch
    TFile hfile(outfile,"RECREATE","simplified ROOT file with single branch structure for training on simulated events");

    // Create a ROOT Tree
    TTree *trainTree = new TTree("Training","ROOT tree with single training branch");
    //trainTree->Branch("triggeredSimEvent", &trigSimEvent, leafList);
    trainTree->Branch("fEnergyGeV", &fEnergyGeV, "fEnergyGeV");
    trainTree->Branch("fLength", &fLength, "fLength");
    trainTree->Branch("fWidth", &fWidth, "fWidth");
    trainTree->Branch("fSize", &fSize, "fSize");
    trainTree->Branch("fDist", &fDist, "fDist");
    trainTree->Branch("fLoss", &fLoss, "fLoss");
    trainTree->Branch("fTimeGradient", &fTimeGradient, "fTimeGradient");
    trainTree->Branch("fPrimaryAzimuthDeg", &fPrimaryAzimuthDeg, "fPrimaryAzimuthDeg");
    // should make 4 branches for each telescope 
    //noise parameter (same as when producing tables)
    
    
    //// set up input to connect simulated events to shower events /////
    
    // open input file for reading only
    VARootIO* pfRootIO = new VARootIO(string(infile.Data()), true); 

    pfRootIO->loadTheRootFile();
    if(!pfRootIO->IsOpen())
        {
            cerr << "Couldn't open " << infile << endl;
	    return 1;
        }

    // parameterised event data holds shower params for each telescope 
    VAParameterisedEventData* paramData = nullptr;
    TTree* paramEventsTree = pfRootIO->loadTheParameterisedEventTree();
    paramEventsTree->SetBranchAddress("P", &paramData);
    //vector<VAParameterisedEventTelData*> vfTels
    //VAHillasData*           pfHillasData;

    // get the simulation tree for true energy, etc. 
    TTree* simTree = pfRootIO->loadTheSimulationEventTree();    
    VASimulationData* simData = nullptr;
    simTree->SetBranchAddress("Sim", &simData);

    // check for trees and branches 
    if(!paramEventsTree)
      {
	  cerr << "Couldn't get a pointer to the shower tree" << endl;
	  return 1;
      }
    if(!simTree)
      {
	  cerr << "Couldn't get a pointer to the simulation tree" << endl;
	  return 1;
      }
    if(!simTree->GetBranch(gSimulatedEventsBranchName.c_str()))
	{
	    cerr << "Couldn't find a branch called " << gSimulatedEventsBranchName << endl;
	    return false;
	}

    
    cout << "simulated events: " << simTree->GetEntries() << endl;
    cout << "parameterised events: " << paramEventsTree->GetEntries() << endl;

    cout << "Building simulation index" << endl;
    paramEventsTree->BuildIndex("fRunNum", "fArrayEventNum");    
    //paramEventsTree->BuildIndex("P.fRunNum", "P.fArrayEventNum");
    //simTree->BuildIndex("fRunNum", "fArrayEventNum");
    cout << "done" << endl;
    

    // now loop over all entries 
    unsigned long long int simEntries = simTree->GetEntries();    
    for(unsigned long long int i = 0; i < simEntries; i++)
      {
	  if(simTree->GetEntry(i) <= 0)
	    {
		cerr << "could not get entry " << i <<endl;
		return 1;
	    }

	  if(i%1000000 == 0)
	  {
	      cout << "event\t" << i <<endl;
//.........这里部分代码省略.........

void MCReader(int eventNo, TString option) {
    int tmpcode;

    cout << "start: " << flush;
    tmpcode = system("date");

    DSReader r("test3.root");
    DS::Root *ds = r.GetEvent(eventNo-3);
    cout << "ds loaded: " << flush;
    tmpcode = system("date");

    DS::MC *mc = ds->GetMC();
    int size_mc = mc->GetMCTrackCount();
    cout << "total tracks: " << size_mc << endl;

    vector<float> xx, yy, zz, tt, ll, wlwl;

    for (int i=0; i < size_mc; i++) {
        DS::MCTrack *track = mc->GetMCTrack(i);
        if (i % 100000 == 0) cout << "processing track " << i << endl;
        if (! (track->GetParticleName() == "opticalphoton") ) continue;
        int size_step = track->GetMCTrackStepCount();
        for (int j=0; j < size_step; j++) {
            DS::MCTrackStep *step = track->GetMCTrackStep(j); // Copy over initial step
            if ( option == "Cerenkov" ||
                    option == "Transportation" ||
                    option == "Attenuation" ||
                    option == "G4FastSimulationManagerProcess"
               ) {
                if (! (step->GetProcess() == option) ) continue;
            }
            else {
                cout << "Sorry, option not recognized" << endl;
            }
            TVector3 vec = step->GetEndpoint();
            // if (vec.Mag()>10000) continue;
            xx.push_back(vec.X());
            yy.push_back(vec.Y());
            zz.push_back(vec.Z());
            tt.push_back(step->GetGlobalTime());
            ll.push_back(step->GetLength());
            wlwl.push_back(1240e-6/step->GetKE()); // nm

            TString name = step->GetProcess();
            if (! ( name == "Cerenkov" ||
                    name == "Transportation" ||
                    name == "Attenuation" ||
                    name == "G4FastSimulationManagerProcess")
               ) {
                cout << step->GetProcess() << endl;
            }

            // if (! (step->GetProcess() == "Cerenkov") ) {

            // }
        }
    }
    cout << "end processing: " << flush;
    tmpcode = system("date");

    TString rootFileName = TString::Format("%i-%s.root", eventNo, option.Data());
    float x, y, z, t, l, wl;
    TTree *tree = new TTree("MCTrack", "MC Track Info");
    tree->Branch("x", &x);
    tree->Branch("y", &y);
    tree->Branch("z", &z);
    tree->Branch("t", &t);
    tree->Branch("l", &l);
    tree->Branch("wl", &wl);
    TFile f(rootFileName.Data(), "recreate");
    int size = xx.size();
    for (int i=0; i<size; i++) {
        x = xx[i];
        y = yy[i];
        z = zz[i];
        t = tt[i];
        wl = wlwl[i];
        l = ll[i];
        tree->Fill();
    }
    tree->Write();
    // TGraph2D g(x.size(), &x[0], &y[0], &z[0]);
    // g.Write();
    f.Close();
    cout << "ROOT file saved to " << rootFileName << endl;

    // TrackNav nav(ds, true);
    // TrackCursor c = nav.Cursor(true);
    // cout << "track loaded: " << flush;  tmpcode = system("date");

    TString jsonFileName = TString::Format("%i-%s.json", eventNo, option.Data());
    ofstream out(jsonFileName.Data());

    out << fixed << setprecision(0);
    out << "{" << endl;
    print_vector(out, xx, "x");
    print_vector(out, yy, "y");
    print_vector(out, zz, "z");

    out << fixed << setprecision(1);
//.........这里部分代码省略.........

int main( int argc, char* argv[] ) {


  std::string runName = "precalib_BGO_pedestal_noSource";
  if( argc>1 ) {
    std::string runName_str(argv[1]);
    runName = runName_str;
  }

  std::string tag = "V00";
  if( argc>2 ) {
    std::string tag_str(argv[2]);
    tag = tag_str;
  }

  TString runName_tstr(runName);
  bool isOnlyRunNumber = !(runName_tstr.BeginsWith("BTF_"));


  TChain* tree = new TChain("recoTree");
  if( isOnlyRunNumber ) {
    std::cout << "-> We believe you are passing the program only the run number!" << std::endl;
    std::cout << "-> So for instance you are passing '246' for run 'BTF_246_20140501-212512_beam'" << std::endl;
    std::cout << "(if this is not the case this means TROUBLE)" << std::endl;
    std::cout << "-> Will look for runs matching run number: " << runName << std::endl;
    tree->Add(Form("analysisTrees_%s/Reco_BTF_%s_*beam.root/recoTree", tag.c_str(), runName.c_str()) );
    if( tree->GetEntries()==0 ) {
      std::cout << "WARNING! Didn't find any events matching run: " << runName << std::endl;
      std::cout << "Exiting" << std::endl;
      exit(1913);
    }
  } else {
    std::string fileName = "analysisTrees_"+tag+"/Reco_" + runName + ".root";
    TFile* file = TFile::Open(fileName.c_str());
    if( file==0 ) {
      std::cout << "ERROR! Din't find file " << fileName << std::endl;
      std::cout << "Exiting." << std::endl;
      exit(11);
    }
    tree = (TChain*)file->Get("recoTree");
  }

  UInt_t evtNumber;
  tree->SetBranchAddress( "evtNumber", &evtNumber );
  UInt_t adcData[40];
  tree->SetBranchAddress( "adcData", adcData );
  UInt_t adcBoard[40];
  tree->SetBranchAddress( "adcBoard", adcBoard );
  UInt_t adcChannel[40];
  tree->SetBranchAddress( "adcChannel", adcChannel );


  unsigned int runNumber;
  int nHodoFibersX;
  int nHodoFibersY;
  int nHodoClustersX;
  int nHodoClustersY;
  float cef3_corr[CEF3_CHANNELS];
  float bgo_corr[BGO_CHANNELS];
  float scintFront;
  float pos_hodoClustX[HODOX_CHANNELS];
  float pos_hodoClustY[HODOY_CHANNELS];
  int nFibres_hodoClustX[HODOX_CHANNELS];
  int nFibres_hodoClustY[HODOY_CHANNELS];
  float xBeam, yBeam;
  bool isSingleEle_scintFront;
  bool cef3_ok;
  bool cef3_corr_ok;
  bool bgo_ok;
  bool bgo_corr_ok;

  tree->SetBranchAddress( "runNumber", &runNumber );

  tree->SetBranchAddress( "scintFront", &scintFront );
  tree->SetBranchAddress( "cef3_corr", cef3_corr );
  tree->SetBranchAddress( "bgo_corr", bgo_corr );

  tree->SetBranchAddress( "nHodoFibersX", &nHodoFibersX );
  tree->SetBranchAddress( "nHodoFibersY", &nHodoFibersY );
  tree->SetBranchAddress( "nHodoClustersX", &nHodoClustersX );
  tree->SetBranchAddress( "pos_hodoClustX", pos_hodoClustX );
  tree->SetBranchAddress( "nFibres_hodoClustX", nFibres_hodoClustX );
  tree->SetBranchAddress( "nHodoClustersY", &nHodoClustersY );
  tree->SetBranchAddress( "pos_hodoClustY", pos_hodoClustY );
  tree->SetBranchAddress( "nFibres_hodoClustY", nFibres_hodoClustY );
  tree->SetBranchAddress( "scintFront", &scintFront );
  tree->SetBranchAddress( "isSingleEle_scintFront", &isSingleEle_scintFront );
  tree->SetBranchAddress( "xBeam", &xBeam );
  tree->SetBranchAddress( "yBeam", &yBeam );

  tree->SetBranchAddress( "cef3_ok", &cef3_ok );
  tree->SetBranchAddress( "cef3_corr_ok", &cef3_corr_ok );
  tree->SetBranchAddress( "bgo_ok", &bgo_ok );
  tree->SetBranchAddress( "bgo_corr_ok", &bgo_corr_ok );

  int nentries = tree->GetEntries();


  if( isOnlyRunNumber ) {
    // modify runname in such a way that it's useful for getBeamPosition and outfile:
//.........这里部分代码省略.........

//for hands on exercise
void reconstruction(){

  // Open the raw data file
  unsigned int    events(-9999);
  unsigned int    run(-9999);
  unsigned int    tdcPhase(9999);

  std::vector<int> *hits_S0 = 0;
  std::vector<int> *hits_S1 = 0;
  std::vector<int> *hits_S2 = 0;
  std::vector<int> *hits_S3 = 0;

  TFile *dutFile = TFile::Open("MyInputFile.root");
  TTree *duttree = (TTree *) dutFile->Get("tree");
  duttree->SetBranchAddress("event", &events);
  duttree->SetBranchAddress("run", &run);
  duttree->SetBranchAddress("dut0_row", &hits_S0);
  duttree->SetBranchAddress("dut1_row", &hits_S1);
  duttree->SetBranchAddress("dut2_row", &hits_S2);
  duttree->SetBranchAddress("dut3_row", &hits_S3);
  duttree->SetBranchAddress("tdcPhase", &tdcPhase);


  // Book all histograms:
  TH1F *hitmap_S0 = new TH1F("hitmap_S0","hitmap_S0", 256, 0. , 256.);
  TH1F *hitmap_S1 = new TH1F("hitmap_S1","hitmap_S1", 256, 0. , 256.);
  TH1F *hitmap_S2 = new TH1F("hitmap_S2","hitmap_S2", 256, 0. , 256.);
  TH1F *hitmap_S3 = new TH1F("hitmap_S3","hitmap_S3", 256, 0. , 256.);

  TH1F *hitmap_S0_masked = new TH1F("hitmap_S0_masked","hitmap_S0_masked", 256, 0. , 256.);
  TH1F *hitmap_S1_masked = new TH1F("hitmap_S1_masked","hitmap_S1_masked", 256, 0. , 256.);
  TH1F *hitmap_S2_masked = new TH1F("hitmap_S2_masked","hitmap_S2_masked", 256, 0. , 256.);
  TH1F *hitmap_S3_masked = new TH1F("hitmap_S3_masked","hitmap_S3_masked", 256, 0. , 256.);

  TH1F *h_clusters_S2 = new TH1F("clusters_S2","clusters_S2", 256, 0. , 256.);
  TH1F *h_clusters_S3 = new TH1F("clusters_S3","clusters_S3", 256, 0. , 256.);


  // Prepare output Reco ROOT Tree:
  std::vector<int> clusterwidth_S0;
  std::vector<double> clusterposition_S0;
  std::vector<int> clusterwidth_S1;
  std::vector<double> clusterposition_S1;
  std::vector<int> clusterwidth_S2;
  std::vector<double> clusterposition_S2;
  std::vector<int> clusterwidth_S3;
  std::vector<double> clusterposition_S3;
  std::vector<double> stubs_M1;
  std::vector<double> stubs_M2;
  std::vector<int> hits_valid_S0;
  std::vector<int> hits_valid_S1;
  std::vector<int> hits_valid_S2;
  std::vector<int> hits_valid_S3;

  TFile *recofile = new TFile("MyOutFile.root","RECREATE");
  TTree *recotree = new TTree("recotree","reconstructed");
  recotree->Branch("events" , &events );
  recotree->Branch("run" , &run );
  recotree->Branch("row_S0", &hits_S0);
  recotree->Branch("row_S1", &hits_S1);
  recotree->Branch("row_S2", &hits_S2);
  recotree->Branch("row_S3", &hits_S3);
  recotree->Branch("row_S0_masked", &hits_valid_S0);
  recotree->Branch("row_S1_masked", &hits_valid_S1);
  recotree->Branch("row_S2_masked", &hits_valid_S2);
  recotree->Branch("row_S3_masked", &hits_valid_S3);
  recotree->Branch("clusters_S0", &clusterposition_S0);
  recotree->Branch("clusters_S1", &clusterposition_S1);
  recotree->Branch("clusters_S2", &clusterposition_S2);
  recotree->Branch("clusters_S3", &clusterposition_S3);
  recotree->Branch("clusterwidth_S0", &clusterwidth_S0);
  recotree->Branch("clusterwidth_S1", &clusterwidth_S1);
  recotree->Branch("clusterwidth_S2", &clusterwidth_S2);
  recotree->Branch("clusterwidth_S3", &clusterwidth_S3);
  recotree->Branch("stubs_module_1",&stubs_M1);
  recotree->Branch("stubs_module_2",&stubs_M2);
  recotree->Branch("tdc", &tdcPhase);


  // The main events loop:
  for (int i = 0; i < duttree->GetEntries(); i++) {
    // Some printout to see where we are:
    if(i%1000 == 0) std::cout << "Processing event " << i << std::endl;
    duttree->GetEntry(i);

    // Masking noisy strips
    // If a strip is "noisy" can be read from the hit map histograms
    // and results in a hardcoded cut on strip numbers.
    for(size_t hit = 0 ; hit < hits_S0->size(); hit++) hitmap_S0->Fill(hits_S0->at(hit));
    // ...

    // Hint: S0 and S1 have 254 strips each, S2, S3 have 256 strips each.

    // Store all strips which are fine in vectors for further analyses
    // hits_valid_S0;


    // Clustering
    // On each sensor, cluster all strips which belong together. 
//.........这里部分代码省略.........

//void MakeFullBigRIPSTree(char *outfile="bigrips.root"){
void MakeFullBigRIPSTree(char *infile, char *outfile="bigrips.root"){

  //  signal(SIGINT,stop_interrupt); // CTRL + C , interrupt

  gSystem->Load("libXMLParser.so");
  gSystem->Load("libanaroot.so");

  TArtStoreManager * sman = TArtStoreManager::Instance();
  TArtEventStore *estore = new TArtEventStore();
  estore->SetInterrupt(&stoploop);
  estore->Open(infile);
  //  estore->Open(2); // bigrips online

  TArtBigRIPSParameters *para = TArtBigRIPSParameters::Instance();
  para->LoadParameter("db/BigRIPSPPAC.xml");
  para->LoadParameter("db/BigRIPSPlastic.xml");
  para->LoadParameter("db/BigRIPSIC.xml");
  para->LoadParameter("db/FocalPlane.xml");
  TArtCalibPID *brcalib = new TArtCalibPID();
  TArtCalibCoin *calcoin= new TArtCalibCoin();
  TArtRecoPID  *brreco = new TArtRecoPID();

  // setting up rips parameters
  TArtRIPS * rips3to5 = brreco->DefineNewRIPS(3,5,"matrix/mat1.mat",7.5715); // f3 - f5
  TArtRIPS * rips5to7 = brreco->DefineNewRIPS(5,7,"matrix/mat2.mat",7.0675); // f5 - f7

  // setting up tof
  TArtTOF * tof3to7  = brreco->DefineNewTOF("F3pl","F7pl",242.513,5);

  // setting up beam id devices
  //  TArtBeam *beam35 = brreco->DefineNewBeam(rips3to5,tof3to7,"F7IC");
  //  TArtBeam *beam57 = brreco->DefineNewBeam(rips5to7,tof3to7,"F7IC");
  //  TArtBeam *beam35ic11 = brreco->DefineNewBeam(rips3to5,tof3to7,"F11IC");
  //TArtBeam *beam57ic11 = brreco->DefineNewBeam(rips5to7,tof3to7,"F11IC");
  TArtBeam *beam37ic7 = brreco->DefineNewBeam(rips3to5,rips5to7,tof3to7,"F11IC");
  //  TArtBeam *beam37ic11 = brreco->DefineNewBeam(rips3to5,rips5to7,tof3to7,"F11IC");

  TFile *fout = new TFile(outfile,"RECREATE");
  TTree *tree = new TTree("tree","tree");

  Int_t maxbit;
  tree->Branch("maxbit",&maxbit,"maxbit/I");
  // define data nodes which are supposed to be dumped to tree 
  TClonesArray * info_array = 
    (TClonesArray *)sman->FindDataContainer("EventInfo");
  tree->Branch(info_array->GetName(),&info_array);
  TClonesArray * ppac_array = 
    (TClonesArray *)sman->FindDataContainer("BigRIPSPPAC");
  tree->Branch(ppac_array->GetName(),&ppac_array);
  TClonesArray * pla_array = 
    (TClonesArray *)sman->FindDataContainer("BigRIPSPlastic");
  tree->Branch(pla_array->GetName(),&pla_array);
  TClonesArray * ic_array = 
    (TClonesArray *)sman->FindDataContainer("BigRIPSIC");
  tree->Branch(ic_array->GetName(),&ic_array);
  TClonesArray * fpl_array = 
    (TClonesArray *)sman->FindDataContainer("BigRIPSFocalPlane");
  tree->Branch(fpl_array->GetName(),&fpl_array);
  TClonesArray * tof_array = 
    (TClonesArray *)sman->FindDataContainer("BigRIPSTOF");
  tree->Branch(tof_array->GetName(),&tof_array);
  TClonesArray * rips_array = 
    (TClonesArray *)sman->FindDataContainer("BigRIPSRIPS");
  tree->Branch(rips_array->GetName(),&rips_array);
  TClonesArray * beam_array = 
    (TClonesArray *)sman->FindDataContainer("BigRIPSBeam");
  tree->Branch(beam_array->GetName(),&beam_array);

  TArtBeam *beam = (TArtBeam*)beam_array->At(0);

  int neve = 0;
  while(estore->GetNextEvent()){
    maxbit = 0;
    if(neve%1000==0)
      cout << "event: " << neve << endl;

    calcoin->ClearData();
    calcoin->LoadData();
    brcalib->ClearData();
    brcalib->ReconstructData();
    brreco->ClearData();
    brreco->ReconstructData();

    //    beam->SetAoQ(beam->GetAoQ()-0.042);

    Int_t tbit = ((TArtEventInfo*)info_array->At(0))->GetTriggerBit();
    for(int i=1;;i++){
      if(tbit%2 == 1)
	maxbit = i;
      else break;

      tbit = tbit/2;
    }
    tree->Fill();

    neve ++;
  }

  fout->Write();
//.........这里部分代码省略.........

//Make pull plots from the output of a toy study
void ResultFormatter::WriteFlatNtuple( const string FileName, const FitResultVector* ToyResult, const vector<string> inputXML, const vector<string> runtimeArgs, const string XMLForProjections, const string XMLForToys )
{
	TFile * rootFile = TFile::Open( FileName.c_str(), "RECREATE" );
	rootFile->SetCompressionLevel( 9 );

	//cout << "Storing Fit Result Vector" << endl;

	//	Important!
	//	The output from this is typically run through the RapidPlot file
	//	For sake of backwards compatibility the RapidPlot tool makes use of the ability to look for the 'first' TTree in a ROOT file for some of it's internal logic
	//	KEEP THIS TREE AS THE FIRST TREE CREATED AND WRITTEN TO THE FILE TO BE ABLE TO KEEP USING THIS TOOL!!!
	TTree* outputTree = new TTree( "RapidFitResult", "RapidFitResult" );

	ResultParameterSet* resultSet = ToyResult->GetFitResult( 0 )->GetResultParameterSet();

	vector<string> allNames = resultSet->GetAllNames();

	for( unsigned int param_i=0; param_i< allNames.size(); ++param_i )
	{
		string thisParamName( allNames[param_i] );

		vector<double> ParameterValues, ParameterErrors, ParameterOriginalValues;
		vector<double> ParameterPulls, ParameterStepSizes;
		vector<double> ParameterErrorsHigh, ParameterErrorsLow;
		vector<double> ParameterMinimums, ParameterMaximums;
		vector<int> ParameterScanStatus, ParameterFixedStatus;
		for( unsigned int resultNum=0; resultNum< (unsigned)ToyResult->NumberResults(); ++resultNum )
		{
			ResultParameter* thisParam = ToyResult->GetFitResult( (unsigned)resultNum )->GetResultParameterSet()->GetResultParameter( thisParamName );

			ParameterValues.push_back( thisParam->GetValue() );
			ParameterErrors.push_back( thisParam->GetError() );
			ParameterPulls.push_back( thisParam->GetPull() );
			ParameterMinimums.push_back( thisParam->GetMinimum() );
			ParameterMaximums.push_back( thisParam->GetMaximum() );
			ParameterOriginalValues.push_back( thisParam->GetOriginalValue() );

			if( thisParam->GetType() == "Fixed" )
			{
				ParameterFixedStatus.push_back( 1 );
			}
			else
			{
				ParameterFixedStatus.push_back( 0 );
			}

			//if( thisParam->GetScanStatus() ) cout << "Scanned: " << string(thisParamName) << endl;

			if( thisParam->GetScanStatus() )
			{
				ParameterScanStatus.push_back( 1 );
			}
			else
			{
				ParameterScanStatus.push_back( 0 );
			}

			if( thisParam->GetAssym() )
			{
				ParameterErrorsHigh.push_back( thisParam->GetErrHi() );
				ParameterErrorsLow.push_back( thisParam->GetErrLow() );
			} else {
				ParameterErrorsHigh.push_back( 0. );
				ParameterErrorsLow.push_back( 0. );
			}
			ParameterStepSizes.push_back( thisParam->GetStepSize() );
		}
		string BranchName=allNames[param_i];
		ResultFormatter::AddBranch( outputTree, BranchName+"_value", ParameterValues );
		bool fixed_param = ToyResult->GetFitResult(0)->GetResultParameterSet()->GetResultParameter(thisParamName)->GetType() == "Fixed";
		bool scanned_param = ToyResult->GetFitResult(0)->GetResultParameterSet()->GetResultParameter(thisParamName)->GetScanStatus();
		if( (!fixed_param) || (scanned_param) )
		{
			ResultFormatter::AddBranch( outputTree, BranchName+"_error", ParameterErrors );
			ResultFormatter::AddBranch( outputTree, BranchName+"_gen", ParameterOriginalValues );
			ResultFormatter::AddBranch( outputTree, BranchName+"_pull", ParameterPulls );
			ResultFormatter::AddBranch( outputTree, BranchName+"_max", ParameterMaximums );
			ResultFormatter::AddBranch( outputTree, BranchName+"_min", ParameterMinimums );
			ResultFormatter::AddBranch( outputTree, BranchName+"_step", ParameterStepSizes );
			ResultFormatter::AddBranch( outputTree, BranchName+"_errHi", ParameterErrorsHigh );
			ResultFormatter::AddBranch( outputTree, BranchName+"_errLo", ParameterErrorsLow );
		}
		ResultFormatter::AddBranch( outputTree, BranchName+"_scan", ParameterScanStatus );
		ResultFormatter::AddBranch( outputTree, BranchName+"_fix", ParameterFixedStatus );
	}

	//cout << "Stored Parameters" << endl;

	ResultFormatter::AddBranch( outputTree, "Fit_RealTime", ToyResult->GetAllRealTimes() );
	ResultFormatter::AddBranch( outputTree, "Fit_CPUTime", ToyResult->GetAllCPUTimes() );
	ResultFormatter::AddBranch( outputTree, "Fit_GLTime", ToyResult->GetAllGLTimes() );

	vector<int> fitStatus;
	vector<double> NLL_Values, RealTimes, CPUTimes;
	for( unsigned int i=0; i< (unsigned) ToyResult->NumberResults(); ++i )
	{
		fitStatus.push_back( ToyResult->GetFitResult( (int)i )->GetFitStatus() );
		NLL_Values.push_back(  ToyResult->GetFitResult( (int)i )->GetMinimumValue() );
	}
//.........这里部分代码省略.........

int main(int argc, char **argv)
{
	struct DanssEventStruct4		DanssEvent;
	struct DanssCmStruct			DanssCm;
	struct DanssEventStruct4		SavedEvent;
	struct DanssInfoStruct4			DanssInfo;
	struct DanssInfoStruct			SumInfo;

	TChain *EventChain;
	TChain *InfoChain;
	TTree *tOut;
	TTree *InfoOut;
	TFile *fOut;
	FILE *fList;
	char str[1024];
	long long iEvt, nEvt;
	long long lastgTime, lastVeto;
	int CmCnt;
	int i;
	char *ptr;
	int nCnt;
	
	if (argc < 3) {
		printf("Usage: ./cmbuilder list_file.txt output_file.root\n");
		printf("Will process files in the list_file and create root-file\n");
		return 10;
	}

	fOut = new TFile(argv[2], "RECREATE");
	if (!fOut->IsOpen()) {
		printf("Can not open the output file %s: %m\n", argv[2]);
		return -10;
	}

	tOut = new TTree("DanssCm", "Time Correlated events");
	tOut->Branch("Cm", &DanssCm,
		"number[10]/L:"		// event numbers in the file
		"unixTime/I:"		// linux time, seconds
		"N/I:"			// number of neutrons + 1
		"SiPmCleanEnergy[10]/F:"	// Full Clean energy SiPm
		"PmtCleanEnergy[10]/F:"	// Full Clean energy Pmt
//		"neutron" parameters
		"Hits[10]/I:"		// SiPm clean hits
//		"NeutronHits[10]/I:"	// number of hits considered as neutron capture gammas
		"NeutronEnergy[10]/F:"	// Energy sum of above (SiPM)
		"NeutronX[10][3]/F:"	// center of gammas position
		"PositronX[10][3]/F:"	// maximum hit clusters
//		"NeutronGammaEnergy[10][5]/F:"	// sorted list of the 5 most energetic gammas
//		"NeutronGammaDistance[10][5]/F:"	// distances for the gammas above to the "neutron" center
		"PositronEnergy[10]/F:"	// maximum hit clusters
		"NeutronRadius[10]/F:"	// average distance between hits and the center
//		Pair parameters
		"gtDiff[10]/F:"		// time difference in us (from 125 MHz clock)
		"Distance[10]/F:"		// distance between neutron and positron, cm
		"DistanceZ[10]/F"		// in Z, cm
	);
	
	InfoOut = new TTree("SumInfo", "Summary information");
	InfoOut->Branch("Info", &SumInfo,  
		"gTime/L:"		// running time in terms of 125 MHz
		"startTime/I:"		// linux start time, seconds
		"stopTime/I:"		// linux stop time, seconds
		"events/I"		// number of events
	);
	memset(&SumInfo, 0, sizeof(struct DanssInfoStruct4));
	memset(&DanssCm, 0, sizeof(struct DanssCmStruct));

	EventChain = new TChain("DanssEvent");
	EventChain->SetBranchAddress("Data", &DanssEvent);
	InfoChain = new TChain("DanssInfo");
	InfoChain->SetBranchAddress("Info", &DanssInfo);

	fList = fopen(argv[1], "rt");
	if (!fList) {
		printf("Can not open list of files %s: %m\n", argv[1]);
		goto fin;
	}
	
	for(;;) {
		if (!fgets(str, sizeof(str), fList)) break;
		ptr = strchr(str, '\n');
		if (ptr) *ptr = '\0';
		EventChain->Add(str);
		InfoChain->Add(str);
	}
	fclose(fList);

	nEvt = EventChain->GetEntries();
	CmCnt = 0;
	nCnt = 0;
	lastVeto = lastgTime = -GLOBALFREQ;
	for (iEvt =0; iEvt < nEvt; iEvt++) {
		EventChain->GetEntry(iEvt);
		if (IsVeto(&DanssEvent)) {
			lastVeto = DanssEvent.globalTime;
			lastgTime = -GLOBALFREQ;
			nCnt = 0;
		}
		if (DanssEvent.globalTime - lastVeto < VETOBLK * GFREQ2US) {
			continue;	// Veto is active
//.........这里部分代码省略.........

void RecoSAMURAI(char * finname="calib.root", 
		 char * foutname="reco.root",
		 Int_t nevents = 99999999,
		 Int_t starteventn = 0 ){

  Double_t samurai_center_brho = 7.2751; // Tm
  Double_t hodpos_offset = 246.44; // mm
  Double_t hodtof_offset = 146.5; // ns
  Double_t dist_pla2target = 2824.; // mm

  TFile *fin = new TFile(finname,"READ");
  gROOT->cd();
  TTree *itree = (TTree*)fin->Get("sm");
  TArtStoreManager * sman = TArtStoreManager::Instance();
  sman->LoadDataFromTree(itree);

  TArtRecoFragment * recofrag = new TArtRecoFragment;
  recofrag->SetCenterBrho(samurai_center_brho);
  recofrag->SetTOFOffSet(hodtof_offset);
  recofrag->SetPla2TargetLength(dist_pla2target);

  std::cout << "getting data containers" << std::endl;
  TClonesArray *fragments = (TClonesArray *)sman->FindDataContainer("SAMURAIFragment");

  Double_t fragdelta, fragbrho, fragtof, fragbeta, fragaoq;
  Double_t trchi2, trpos[3], trl2hod;
  Int_t trstatus, trnhit;

  TFile *fout = new TFile(foutname,"RECREATE");
  TTree *otree = itree->CloneTree(0);
  /*
  TObjArray *branchArray = itree->GetListOfBranches();
  for (Int_t i=0; i<(size_t)(branchArray->GetEntriesFast()); i++) {
    TBranch *branch = (TBranch *) ((*branchArray)[i]);
    otree->Branch(branch->GetName(),&branch);
  }
  */

  otree->Branch("fragdelta",&fragdelta,"fragdelta/D");
  otree->Branch("fragbrho",&fragbrho,"fragbrho/D");
  otree->Branch("fragtof",&fragtof,"fragtof/D");
  otree->Branch("fragbeta",&fragbeta,"fragbeta/D");
  otree->Branch("fragaoq",&fragaoq,"fragaoq/D");
  otree->Branch("trstatus",&trstatus,"trstatus/I");
  otree->Branch("trchi2",&trchi2,"trchi2/D");
  otree->Branch("trnhit",&trnhit,"trnhit/I");
  otree->Branch("trpos",trpos,"trpos[3]/D");
  otree->Branch("trl2hod",&trl2hod,"trl2hod/D");

  int nentries = itree->GetEntries();
  int neve = starteventn;
  std::cout << "entries in input tree: " << nentries << std::endl;
  std::cout << "start to scan from event" << neve << std::endl;

  while(neve<nevents+starteventn && neve < nentries){

    itree->GetEntry(neve);
    recofrag->ClearData();
    recofrag->ReconstructData();

    if(fragments){
      if(fragments->GetEntriesFast()>0){
	TArtFragment *f= (TArtFragment*)fragments->At(0);
	fragdelta = f->GetDelta();
	fragbrho = f->GetBrho();
	//	std::cout << fragbrho << std::endl;
	fragtof = f->GetTOF();
	fragbeta = f->GetBeta();
	fragaoq = f->GetAoQ();
	trchi2 = f->GetChi2();
	trpos[0] = f->GetPrimaryPosition().X();
	trpos[1] = f->GetPrimaryPosition().Y();
	trpos[2] = f->GetPrimaryPosition().Z();
	trl2hod = f->GetPathLengthToHOD();
	trstatus = f->GetRKtraceStatus() ? 1 : 0;
	trnhit = f->GetNHit();
      }
      else{
	fragdelta = -9999;
	fragbrho = -9999;
	fragtof = -9999;
	fragbeta = -9999;
	fragaoq = -9999;
	trchi2 = -9999;
	for(int i=0;i<3;i++) trpos[i] = -9999;
	trl2hod = -9999;
	trstatus = 0;
	trnhit = -1;
      }
    }
    else{
      fragdelta = -9999;
      fragbrho = -9999;
      fragtof = -9999;
      fragbeta = -9999;
      fragaoq = -9999;
      trchi2 = -9999;
      for(int i=0;i<3;i++) trpos[i] = -9999;
      trl2hod = -9999;
      trstatus = 0;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	// Book output histograms
   	TH1F* hCHFb0_OpenSelection= new TH1F		("hCHFb0_OpenSelection", "charged Hadron Energy Fraction b1", 40, 0.0, 1.2);
   	TH1F* hCHFb1_OpenSelection= new TH1F		("hCHFb1_OpenSelection", "charged Hadron Energy Fraction b2", 40, 0.0, 1.2);
 	TH1F* hPtjj_OpenSelection= new TH1F		("hPtjj_OpenSelection","Pt of two b jets with highest CSV ", 50, 0.0, 400);
 	TH1F* hPtmumu_OpenSelection= new TH1F	("hPtmumu_OpenSelection","Pt of two muons with highest pt ", 50, 0.0, 400);
	TH1F* hPtbalZH_OpenSelection= new TH1F	("hPtbalZH_OpenSelection", "Pt balance of Z and H", 40, -80, 80);
	TH1F* hPtmu0_OpenSelection= new TH1F		("hPtmu0_OpenSelection","Pt of muon with highest pt ", 30, 0.0, 300);
	TH1F* hPtmu1_OpenSelection= new TH1F		("hPtmu1_OpenSelection","Pt of muon with second highest pt ", 30, 0.0, 300);
	TH1F* hPFRelIsomu0_OpenSelection= new TH1F		("hPFRelIsomu0_OpenSelection", "PF Rel Iso of muon with highest Pt", 40, 0., 0.2);
	TH1F* hPFRelIsomu1_OpenSelection= new TH1F		("hPFRelIsomu1_OpenSelection", "PF Rel Iso of muon with second highest Pt", 40, 0., 0.2);
	TH1F* hCSV0_OpenSelection= new TH1F		("hCSV0_OpenSelection","Jet with highest CSV ",			40, 0, 1.5);
	TH1F* hCSV1_OpenSelection= new TH1F		("hCSV1_OpenSelection","Jet with second highest CSV ",		40, 0, 1.5);
	TH1F* hdphiVH_OpenSelection= new TH1F	("hdphiVH_OpenSelection","Delta phi between Z and Higgs ", 50, -0.1, 4.5);
	TH1F* hdetaJJ_OpenSelection= new TH1F	("hdetaJJ_OpenSelection","Delta eta between two jets ", 60, -4, 4);
	TH1F* hNjets_OpenSelection= new TH1F	("hNjets_OpenSelection", "Number of Jets",		13, -2.5, 10.5);
	TH1F* hMjj_OpenSelection	= new TH1F	("hMjj_OpenSelection",  "Invariant Mass of two Jets ",		50, 0, 300);
	TH1F* hMmumu_OpenSelection	= new TH1F	("hMmumu_OpenSelection",  "Invariant Mass of two muons ",	75, 0, 200);
    TH1F* hRMSeta_OpenSelection= new TH1F	("hRMSeta_OpenSelection", "RMS Eta",		30, 0, 3);
    TH1F* hStaDeveta_OpenSelection= new TH1F	("hStaDeveta_OpenSelection", "Standard Deviation Eta",		30, 0, 3);
	TH1F* hUnweightedEta_OpenSelection= new TH1F	("hUnweightedEta_OpenSelection",  "Unweighted Eta ",		50, 0, 15);		
	TH1F* hdphiJJ_vect_OpenSelection= new TH1F	("hdphiJJ_vect_OpenSelection", "Delta phi between two jets",  30, -3.5, 4);
	TH1F* hCircularity_OpenSelection= new TH1F("hCircularity_OpenSelection","EventShapeVariables circularity", 30, 0.0, 1.2);
	TH1F* hHt_OpenSelection= new TH1F("hHt_OpenSelection","scalar sum of pt of four particles", 50, 0.0, 500);
    TH1F* hCentrality_OpenSelection= new TH1F	("hCentrality_OpenSelection", "Centrality", 40, 0.0, 0.8);
    TH1F* hEventPt_OpenSelection= new TH1F	("hEventPt_OpenSelection", "Pt of HV system", 50, 0.0, 100);
    TH1F* hAngle_OpenSelection= new TH1F	("hAngle_OpenSelection", "Angle between H and Z", 45, 0, 4.5);
    TH1F* hSphericity_OpenSelection= new TH1F	("hSphericity_OpenSelection", "EventShapeVariables sphericity", 50, 0.0, 1);
    TH1F* hAplanarity_OpenSelection= new TH1F	("hAplanarity_OpenSelection", "EventShapeVariables Aplanarity", 50, -0.1, .4);
    TH1F* hIsotropy_OpenSelection= new TH1F	("hIsotropy_OpenSelection",  "EventShapeVariables isotropy", 30, 0.0, 1.3);
    TH2F* hDphiDetajj_OpenSelection= new TH2F	("hDphiDetajj_OpenSelection", "#Delta#phi vs #Delta#eta JJ", 25, -5, 5, 25, -5, 5);
	
	TTree *treeWithBDT = new TTree("treeWithBDT","Tree wiht BDT output");
	treeWithBDT->SetDirectory(0);
	treeWithBDT->Branch("nJets",&nJets, "nJets/I");
	treeWithBDT->Branch("Naj",&Naj, "Naj/F");
	treeWithBDT->Branch("eventFlavor",&eventFlavor, "eventFlavor/I");
	treeWithBDT->Branch("CSV0",&CSV0, "CSV0/F");
	treeWithBDT->Branch("CSV1",&CSV1, "CSV1/F");
	treeWithBDT->Branch("Zmass",&Zmass, "Zmass/F");
	treeWithBDT->Branch("Hmass",&Hmass, "Hmass/F");
	treeWithBDT->Branch("DeltaPhiHV",&DeltaPhiHV, "DeltaPhiHV/F");
	treeWithBDT->Branch("Hpt",&Hpt, "Hpt/F");
	treeWithBDT->Branch("Zpt",&Zpt, "Zpt/F");
	treeWithBDT->Branch("mu0pt",&mu0pt, "mu0pt/F");
	treeWithBDT->Branch("Ht",&Ht, "Ht/F");
	treeWithBDT->Branch("EtaStandDev",&EtaStandDev, "EtaStandDev/F");
	treeWithBDT->Branch("UnweightedEta",&UnweightedEta, "UnweightedEta/F");
	treeWithBDT->Branch("EvntShpCircularity",&EvntShpCircularity, "EvntShpCircularity/F");
	treeWithBDT->Branch("alpha_j",&alpha_j, "alpha_j/F");
	treeWithBDT->Branch("qtb1",&qtb1, "qtb1/F");
	treeWithBDT->Branch("nSV",&nSV, "nSV/F");
	treeWithBDT->Branch("Trigweight",&Trigweight, "Trigweight/F");
	treeWithBDT->Branch("B2011PUweight",&B2011PUweight, "B2011PUweight/F");
	treeWithBDT->Branch("btag2CSF",&btag2CSF, "btag2CSF/F");
	treeWithBDT->Branch("DetaJJ",&DetaJJ, "DetaJJ/F");
	treeWithBDT->Branch("jetCHF0",&jetCHF0, "jetCHF0/F");
	treeWithBDT->Branch("jetCHF1",&jetCHF1, "jetCHF1/F");
	treeWithBDT->Branch("jetPhi0",&jetPhi0, "jetPhi0/F");
	treeWithBDT->Branch("jetPhi1",&jetPhi1, "jetPhi1/F");
	treeWithBDT->Branch("jetEta0",&jetEta0, "jetEta0/F");
	treeWithBDT->Branch("jetEta1",&jetEta1, "jetEta1/F");
	treeWithBDT->Branch("mu1pt",&mu1pt, "mu1pt/F");
	treeWithBDT->Branch("muonPFiso0",&muonPFiso0, "muonPFiso0/F");
	treeWithBDT->Branch("muonPFiso1",&muonPFiso1, "muonPFiso1/F");
	treeWithBDT->Branch("DphiJJ",&DphiJJ, "DphiJJ/F");
	treeWithBDT->Branch("RMS_eta",&RMS_eta, "RMS_eta/F");

//.........这里部分代码省略.........
  TClonesArray *genPartArr     = new TClonesArray("baconhep::TGenParticle");
  TClonesArray *electronArr    = new TClonesArray("baconhep::TElectron");
  TClonesArray *scArr          = new TClonesArray("baconhep::TPhoton");
  TClonesArray *vertexArr      = new TClonesArray("baconhep::TVertex");

  TFile *infile=0;
  TTree *eventTree=0;
  
  //
  // loop over samples
  //  
  for(UInt_t isam=0; isam<samplev.size(); isam++) {
    
    // Assume data sample is first sample in .conf file
    // If sample is empty (i.e. contains no ntuple files), skip to next sample
    Bool_t isData=kFALSE;
    if(isam==0 && !hasData) continue;
    else if (isam==0) isData=kTRUE;
    
    // Assume signal sample is given name "zee" - flag to store GEN Z kinematics
    Bool_t isSignal = (snamev[isam].CompareTo("zee",TString::kIgnoreCase)==0);  
    // flag to reject Z->ee events when selecting at wrong-flavor background events
    Bool_t isWrongFlavor = (snamev[isam].CompareTo("zxx",TString::kIgnoreCase)==0);  
    
    CSample* samp = samplev[isam];
  
    //
    // Set up output ntuple
    //
    TString outfilename = ntupDir + TString("/") + snamev[isam] + TString("_select.root");
    if(isam!=0 && !doScaleCorr) outfilename = ntupDir + TString("/") + snamev[isam] + TString("_select.raw.root");
    TFile *outFile = new TFile(outfilename,"RECREATE"); 
    TTree *outTree = new TTree("Events","Events");
    outTree->Branch("runNum",     &runNum,     "runNum/i");      // event run number
    outTree->Branch("lumiSec",    &lumiSec,    "lumiSec/i");     // event lumi section
    outTree->Branch("evtNum",     &evtNum,     "evtNum/i");      // event number
    outTree->Branch("matchGen",   &matchGen,   "matchGen/i");    // event has both leptons matched to MC Z->ll
    outTree->Branch("category",   &category,   "category/i");    // dilepton category
    outTree->Branch("id_1",       &id_1,       "id_1/i");        // PDF info -- parton ID for parton 1
    outTree->Branch("id_2",       &id_2,       "id_2/i");        // PDF info -- parton ID for parton 2
    outTree->Branch("x_1",        &x_1,        "x_1/d");         // PDF info -- x for parton 1
    outTree->Branch("x_2",        &x_2,        "x_2/d");         // PDF info -- x for parton 2
    outTree->Branch("xPDF_1",     &xPDF_1,     "xPDF_1/d");      // PDF info -- x*F for parton 1
    outTree->Branch("xPDF_2",     &xPDF_2,     "xPDF_2/d");      // PDF info -- x*F for parton 2
    outTree->Branch("scalePDF",   &scalePDF,   "scalePDF/d");    // PDF info -- energy scale of parton interaction
    outTree->Branch("weightPDF",  &weightPDF,  "weightPDF/d");   // PDF info -- PDF weight
    outTree->Branch("npv",        &npv,        "npv/i");         // number of primary vertices
    outTree->Branch("npu",        &npu,        "npu/i");         // number of in-time PU events (MC)
    outTree->Branch("genV",      "TLorentzVector",  &genV);      // GEN boson 4-vector
    outTree->Branch("genVPt",     &genVPt,     "genVPt/F");      // GEN boson pT (signal MC)
    outTree->Branch("genVPhi",    &genVPhi,    "genVPhi/F");     // GEN boson phi (signal MC)
    outTree->Branch("genVy",      &genVy,      "genVy/F");       // GEN boson rapidity (signal MC)
    outTree->Branch("genVMass",   &genVMass,   "genVMass/F");    // GEN boson mass (signal MC)
    outTree->Branch("genWeight",   &genWeight,  "genWeight/F");
    outTree->Branch("PUWeight",    &PUWeight,   "PUWeight/F");
    outTree->Branch("scale1fb",   &scale1fb,   "scale1fb/F");    // event weight per 1/fb (MC)
    outTree->Branch("scale1fbUp",    &scale1fbUp,   "scale1fbUp/F");    // event weight per 1/fb (MC)
    outTree->Branch("scale1fbDown",    &scale1fbDown,   "scale1fbDown/F");    // event weight per 1/fb (MC)
    outTree->Branch("met",        &met,        "met/F");         // MET
    outTree->Branch("metPhi",     &metPhi,     "metPhi/F");      // phi(MET)
    outTree->Branch("sumEt",      &sumEt,      "sumEt/F");       // Sum ET
    outTree->Branch("u1",         &u1,         "u1/F");          // parallel component of recoil
    outTree->Branch("u2",         &u2,         "u2/F");          // perpendicular component of recoil
    outTree->Branch("tkMet",      &tkMet,      "tkMet/F");       // MET (track MET)
    outTree->Branch("tkMetPhi",   &tkMetPhi,   "tkMetPhi/F");    // phi(MET) (track MET)
    outTree->Branch("tkSumEt",    &tkSumEt,    "tkSumEt/F");     // Sum ET (track MET)

void AddDetectorEffects(TString infile,TString outfile)
{
	InitExterns();
	const int kMaxTrack = 10;	// maximal number of particles in a single event.
	
	
	int evt = -999;
	
	int nPhoton = 0;
	float photonPt[kMaxTrack];
	float photonEta[kMaxTrack];
	float photonPhi[kMaxTrack];
	
	int nElectron = 0;
	float electronPt[kMaxTrack];
	float electronEta[kMaxTrack];
	float electronPhi[kMaxTrack];
	int electronQ[kMaxTrack];
	int electronTag[kMaxTrack];
	float electronCal02[kMaxTrack];
	float electronCal03[kMaxTrack];
	float electronCal04[kMaxTrack];
	float electronTrk02[kMaxTrack];
	float electronTrk03[kMaxTrack];
	float electronTrk04[kMaxTrack];
	float electronSign[kMaxTrack];
	
	int nMuon = 0;
	float muonPt[kMaxTrack];
	float muonEta[kMaxTrack];
	float muonPhi[kMaxTrack];
	int muonQ[kMaxTrack];
	int muonTag[kMaxTrack];
	float muonCal02[kMaxTrack];
	float muonCal03[kMaxTrack];
	float muonCal04[kMaxTrack];
	float muonTrk02[kMaxTrack];
	float muonTrk03[kMaxTrack];
	float muonTrk04[kMaxTrack];
	float muonSign[kMaxTrack];
	
	int nTau = 0;
	float tauPt[kMaxTrack];
	float tauEta[kMaxTrack];
	float tauPhi[kMaxTrack];
	int tauQ[kMaxTrack];
	int tauTag[kMaxTrack];
	float tauCal02[kMaxTrack];
	float tauCal03[kMaxTrack];
	float tauCal04[kMaxTrack];
	float tauTrk02[kMaxTrack];
	float tauTrk03[kMaxTrack];
	float tauTrk04[kMaxTrack];
	float tauSign[kMaxTrack];
	
	int nJet = 0;
	float jetPt[kMaxTrack];
	float jetEta[kMaxTrack];
	float jetPhi[kMaxTrack];
	int jetQ[kMaxTrack];
	int jetTag[kMaxTrack];
	float jetCal02[kMaxTrack];
	float jetCal03[kMaxTrack];
	float jetCal04[kMaxTrack];
	float jetTrk02[kMaxTrack];
	float jetTrk03[kMaxTrack];
	float jetTrk04[kMaxTrack];
	
	float metPt = -999.;
	float metEta = -999.;
	float metPhi = -999.;

	float sf_mu_reco_eff[kMaxTrack];	// scale factor muon reco efficiency
	float sf_el_reco_eff[kMaxTrack];	// scale factor electron reco efficiency
	
	
	TFile f(outfile,"recreate");	// the file and tree
	TTree *t = new TTree("t","Reconst ntuple");
	
	
	// adding branches:
	// t->Branch(Branch name in the browser, pointer to the variable, variable name)
	
	t->Branch("evt",&evt,"evt/I");

	t->Branch("nPhoton",&nPhoton,"nPhoton/I");
	t->Branch("photonPt",photonPt,"photonPt[nPhoton]/F");
	t->Branch("photonEta",photonEta,"photonEta[nPhoton]/F");
	t->Branch("photonPhi",photonPhi,"photonPhi[nPhoton]/F");

	t->Branch("nElectron",&nElectron,"nElectron/I");
	t->Branch("electronPt",electronPt,"electronPt[nElectron]/F");
	t->Branch("electronEta",electronEta,"electronEta[nElectron]/F");
	t->Branch("electronPhi",electronPhi,"electronPhi[nElectron]/F");
	t->Branch("electronQ",electronQ,"electronQ[nElectron]/I");
	t->Branch("electronTag",electronTag,"electronTag[nElectron]/I");
	t->Branch("electronCal02",electronCal02,"electronCal02[nElectron]/F");
	t->Branch("electronCal03",electronCal03,"electronCal03[nElectron]/F");
	t->Branch("electronCal04",electronCal04,"electronCal04[nElectron]/F");
	t->Branch("electronTrk02",electronTrk02,"electronTrk02[nElectron]/F");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  int firstEvent;
  int initialEvents;
  anaEL.LoadAnalysisInfo("Run", "firstEvent", firstEvent);
  anaEL.LoadAnalysisInfo("Run", "nEvents", nEvents);
  anaEL.LoadAnalysisInfo("Run", "initialEvents", initialEvents);
  string channelName;
  anaEL.LoadAnalysisInfo("MassAnalysis", "ChannelName", channelName);

  anaEL.LoadWeight (sel); // now the parameters for SFBweight are initialized (for b-tag!)

  //******************************************
  //Load Scale Factors for lepton efficiencies
  //******************************************
  sel.LoadElScaleFactors();
  sel.LoadMuScaleFactors();
  sel.InitJESUnc();

  bool nonskimmed = true;

  TopTree::NTEvent * event = 0;

  TFile * TupleFile = new TFile("TupleFileNonSkimmed_TTbarSignalPowheg.root","RECREATE");
  float mttbar_rec, mttbar_gen;
  float diffysquare_rec, diffysquare_gen;
  float diffabseta_rec, diffabseta_gen;
  float weight;
  float pt_ttbar, pt_ttbar_gen;
  float pz_ttbar, pz_ttbar_gen;
  float eta_ttbar_rec, eta_ttbar_gen;
  float y_ttbar_rec, y_ttbar_gen;
  float diffabsy_rec, diffabsy_gen;

  TTree *tree = new TTree("Tuple","Tuple");
  tree->Branch("mttbar_rec",&mttbar_rec,"mttbar_rec/F");
  tree->Branch("mttbar_gen",&mttbar_gen,"mttbar_gen/F");
  tree->Branch("diffysquare_rec",&diffysquare_rec,"diffysquare_rec/F");
  tree->Branch("diffysquare_gen",&diffysquare_gen,"diffysquare_gen/F");
  tree->Branch("diffabseta_rec",&diffabseta_rec,"diffabseta_rec/F");
  tree->Branch("diffabseta_gen",&diffabseta_gen,"diffabseta_gen/F");
  tree->Branch("weight",&weight,"weight/F");
  tree->Branch("pt_ttbar",&pt_ttbar,"pt_ttbar/F");
  tree->Branch("pt_ttbar_gen",&pt_ttbar_gen,"pt_ttbar_gen/F");
  tree->Branch("pz_ttbar",&pz_ttbar,"pz_ttbar/F");
  tree->Branch("pz_ttbar_gen",&pz_ttbar_gen,"pz_ttbar_gen/F");
  tree->Branch("eta_ttbar_rec",&eta_ttbar_rec,"eta_ttbar_rec/F");
  tree->Branch("eta_ttbar_gen",&eta_ttbar_gen,"eta_ttbar_gen/F");
  tree->Branch("y_ttbar_rec",&y_ttbar_rec,"y_ttbar_rec/F");
  tree->Branch("y_ttbar_gen",&y_ttbar_gen,"y_ttbar_gen/F");
  tree->Branch("diffabsy_rec",&diffabsy_rec,"diffabsy_rec/F");
  tree->Branch("diffabsy_gen",&diffabsy_gen,"diffabsy_gen/F");

  //////////////////////
  //LOOP OVER THE DATASETS
  //////////////////////
  if(verbosity>0) {
    cout<<"#########################"<<endl;
    cout<<" Loop over the datasets  "<<endl;
    cout<<"#########################"<<endl;
  }

  float num_mttbar_eff = 0.;
  float den_mttbar_eff = 0.;

  for (unsigned int d = 0; d < datasets.size (); d++) {

    weight_MC = datasets[d].NormFactor()*Luminosity;

//.........这里部分代码省略.........
  Int_t QQTrkPt04;
  Int_t QQTrkPt03;
  Int_t QQTrkPt02;

  float weight;
  float weight2;
  Float_t upsPt;
  Float_t upsEta;
  Float_t upsPhi;
  Float_t upsRapidity;	
  Float_t vProb;
  Float_t muPlusPt;
  Float_t muMinusPt;
  Float_t muPlusEta;
  Float_t muMinusEta;
  Float_t muPlusPhi;
  Float_t muMinusPhi;
  Float_t RmuPlusPhi = 0;
  Float_t RmuMinusPhi = 0;

  // centrliaty extra stuff
  Int_t Npix, NpixelTracks,Ntracks;
  Float_t SumET_HF, SumET_HFplus, SumET_HFminus, SumET_HFplusEta4, SumET_HFminusEta4, SumET_EB, SumET_ET, SumET_EE, SumET_EEplus, SumET_EEminus, SumET_ZDC, SumET_ZDCplus, SumET_ZDCminus;

  // track extra stuff
  const int NTRKMAX=100000000;
  int Reco_trk_size=0;
  TClonesArray *Reco_trk_vtx = 0;
  TClonesArray *Reco_trk_4mom = 0;
  float trkPt[NTRKMAX];
  float trkEta[NTRKMAX];
  float trkPhi[NTRKMAX];
  //---------------------------
  TBranch *b_Centrality; //!
  TBranch *b_eventNb;    //!
  TBranch *b_runNb;      //!
  TBranch *b_HLTriggers; //!
  TBranch *b_Reco_QQ_size; //!
  TBranch *b_Reco_QQ_trig; //!
  TBranch *b_Reco_QQ_type; //!
  TBranch *b_Reco_QQ_sign; //!
  TBranch *b_Reco_QQ_VtxProb;  //!
  TBranch *b_Reco_QQ_dca; // new float, unused in upsilon
  TBranch *b_Reco_QQ_ctau; //!
  TBranch *b_Reco_QQ_ctauErr; //!
  TBranch *b_Reco_QQ_ctauTrue; //!
  TBranch *b_Reco_QQ_4mom;     //!
  TBranch *b_Reco_QQ_vtx;      //!
  TBranch *b_Reco_QQ_mupl_4mom;//!
  TBranch *b_Reco_QQ_mumi_4mom;//!
  TBranch *b_Reco_QQ_NtrkDeltaR03;//!
  TBranch *b_Reco_QQ_NtrkDeltaR04;//!
  TBranch *b_Reco_QQ_NtrkDeltaR05;//!
  TBranch *b_Reco_QQ_NtrkPt04;//!
  TBranch *b_Reco_QQ_NtrkPt03;//!
  TBranch *b_Reco_QQ_NtrkPt02;//!
  //------------------------------
  TBranch *b_Reco_QQ_mupl_nTrkHits;           
  TBranch *b_Reco_QQ_mupl_normChi2_inner;
  TBranch *b_Reco_QQ_mupl_normChi2_global;
  TBranch *b_Reco_QQ_mupl_dxy;
  TBranch *b_Reco_QQ_mupl_dxyErr;
  TBranch *b_Reco_QQ_mupl_dz;
  TBranch *b_Reco_QQ_mupl_dzErr;
  TBranch *b_Reco_QQ_mupl_TrkMuArb;
  TBranch *b_Reco_QQ_mupl_TMOneStaTight; // new bool, unused in upsilon(?)

void splitter() { 
  
  //read input file
  TFile input_file("daq_converted.root","READ");

  //read input trees
  TTree *volumetree = (TTree*)input_file.Get("TestVolumeTree");
  TTree *input_tree0 = (TTree*)input_file.Get("Mod0Tree");
  TTree *input_tree1 = (TTree*)input_file.Get("Mod1Tree");
  TTree *input_tree2 = (TTree*)input_file.Get("Mod2Tree");
  TTree *input_tree3 = (TTree*)input_file.Get("Mod3Tree");

  int nentries0 = input_tree0->GetEntries();
  std::clog << " input has " << nentries0 << " entries" << std::endl;
  
  //associate branches to input trees
  double gap;
  volumetree->SetBranchAddress("Gap",&gap);
  
  std::vector<int>HitPinsTop0;
  input_tree0->SetBranchAddress("HitPinsTop0",&(&HitPinsTop0));
  std::vector<int>HitPinsBot0;
  input_tree0->SetBranchAddress("HitPinsBot0",&(&HitPinsBot0));

  std::vector<int>HitPinsTop1;
  input_tree1->SetBranchAddress("HitPinsTop1",&(&HitPinsTop1));
  std::vector<int>HitPinsBot1;
  input_tree1->SetBranchAddress("HitPinsBot1",&(&HitPinsBot1));

  std::vector<int>HitPinsTop2;
  input_tree2->SetBranchAddress("HitPinsTop2",&(&HitPinsTop2));
  std::vector<int>HitPinsBot2;
  input_tree2->SetBranchAddress("HitPinsBot2",&(&HitPinsBot2));
  
  std:vector<int>HitPinsTop3;
  input_tree3->SetBranchAddress("HitPinsTop3",&(&HitPinsTop3));
  std:vector<int>HitPinsBot3;
  input_tree3->SetBranchAddress("HitPinsBot3",&(&HitPinsBot3));

  //loop over input trees
  int Nentries = 10000;
  int Nfiles = floor(nentries0/Nentries + 0.5);
  
  for(int iFile = 0; iFile<=Nfiles; iFile++){ 
      // make output file
      TFile output_file(Form("daq_split_%i.root",iFile),"RECREATE");
      
      // make output trees
      TTree *modtree0 = new TTree("Mod0Tree","Mod0Tree");
      TTree *modtree1 = new TTree("Mod1Tree","Mod1Tree");
      TTree *modtree2 = new TTree("Mod2Tree","Mod2Tree");
      TTree *modtree3 = new TTree("Mod3Tree","Mod3Tree");
      
      TTree *testvoltree = new TTree("TestVolumeTree","TestVolumeTree");
      double gapn = 411;
      testvoltree->Branch("Gap",&gapn);
      testvoltree->Fill();
      testvoltree->Write();
      
      //create branches for output trees
      
      std::vector<int>HitPinsTop0n;
      modtree0->Branch("HitPinsTop0",&(&HitPinsTop0));
      std::vector<int>HitPinsBot0n;
      modtree0->Branch("HitPinsBot0",&(&HitPinsBot0)); 

      std::vector<int>HitPinsTop1n;
      modtree1->Branch("HitPinsTop1",&(&HitPinsTop1));
      std::vector<int>HitPinsBot1n;
      modtree1->Branch("HitPinsBot1",&(&HitPinsBot1));

      std::vector<int>HitPinsTop2n;
      modtree2->Branch("HitPinsTop2",&(&HitPinsTop2));
      std::vector<int>HitPinsBot2n;
      modtree2->Branch("HitPinsBot2",&(&HitPinsBot2));

      std::vector<int>HitPinsTop3n;
      modtree3->Branch("HitPinsTop3",&(&HitPinsTop3));
      std::vector<int>HitPinsBot3n;
      modtree3->Branch("HitPinsBot3",&(&HitPinsBot3));
     
      //loop for splitting entries
      for(int i = 0; i<Nentries; i++){ 

	int ientry = Nentries*iFile + i ;
	
        if(ientry<nentries0){

	  input_tree0->GetEntry(ientry);
	  modtree0->Fill();

	  input_tree1->GetEntry(ientry);
	  modtree1->Fill();

	  input_tree2->GetEntry(ientry);
	  modtree2->Fill();

	  input_tree3->GetEntry(ientry);
	  modtree3->Fill();
	  	  
//.........这里部分代码省略.........