C++ TBranch::Fill方法代码示例

void updateweight(TString filename)
{
  auto f = new TFile(filename,"update");

  auto nt = (TTree *)f->Get("nt");

  float prew, weight;
  TBranch *bw;

  bw =  nt->Branch("weight",&weight);
  nt->SetBranchAddress("prew",&prew);
  
  int n = nt->GetEntries();
  int onep = n/100;
  for (int i=0;i<n;i++) {
    if (i%onep==0) cout<<i/onep<<endl;
    nt->GetEntry(i);
    weight = prew;
    bw->Fill();
  }

  nt->Write();
  f->Close();


}

void ResultFormatter::AddBranch( TTree* inputTree, const string& BranchName, const vector<double>& DoubleData )
{
	if( inputTree->GetEntries() != 0 )
	{
		if( inputTree->GetEntries() != (int) DoubleData.size() )
		{
			cerr << "CANNOT ADD DOUBLE BRANCH: " << BranchName << " TO: " << inputTree->GetName() << endl;
			return;
		}
	}

	TString branchName( BranchName );
	TString branchLabel=branchName; branchLabel.Append("/D");

	double thisValue=0.;

	TBranch* newBranch = inputTree->Branch( branchName, &thisValue, branchLabel );

	inputTree->SetEntries( (int) DoubleData.size() );

	for( unsigned int i=0; i< DoubleData.size(); ++i )
	{
		thisValue=DoubleData[i];
		newBranch->Fill();
	}
	return;
}

int reader_wrapper::GetEntry(Long64_t e) {
  /// don't care about spectators here
  for (auto b: m_branches) {
    b->GetEntry(e);
  }
  for (auto& v : m_variables) {
    v.value = v.ttreeformula->EvalInstance();
  }
  m_response = m_reader->EvaluateMVA(m_methodName.Data());
  m_responseBranch->Fill();
  return 0;
}

void addBranchToTreesInFiles() {
  //************************************************************
  //                      Variables                           //
  vector<TString> fileName;
  fileName.push_back( "PhotonJetPt15_realPhotons.root" );
  fileName.push_back( "PhotonJetPt30_realPhotons.root" );

  // The following 4 variables set the scale
  float invLuminosityToScaleTo = 200; // in pb-1

  vector<float> crossSection;
  crossSection.push_back( 2.9E5 );
  crossSection.push_back( 3.2E4 );

  vector<float> filterEffeciency;
  filterEffeciency.push_back( 1.0 );
  filterEffeciency.push_back( 1.0 );

  vector<float> eventsAnalyzied;
  eventsAnalyzied.push_back( 950E3 );
  eventsAnalyzied.push_back( 670E3 );
  // END of setting scale

  TString treeName = "TreePhotonMatched";

  //                  END of Variables                        //
  //************************************************************


  // Loop over all the Files
  for (int i=0; i < fileName.size(); i++) {
    TFile* currentFile = new TFile(fileName[i],"update");

    // this is the variable to be added
    Float_t scale=crossSection[i]*invLuminosityToScaleTo*filterEffeciency[i]/eventsAnalyzied[i];

    cout << "Opened " << fileName[i] << ", adding new branch with value=" << scale;

    TTree *tree = (TTree*)currentFile->Get(treeName);
    TBranch *newBranch = tree->Branch("weight", &scale,"weight/F");

    // Loop over all the entries, and add the new branch
    Int_t numEntries = (Int_t)tree->GetEntries();
    for (Int_t j=0; j<numEntries; j++) {
      newBranch->Fill();
    }
    tree->Write("",TObject::kOverwrite); // save new version only
    currentFile->Close();
    cout << "...closed file." << endl;
  }
}

void addDiMuMassPrimeBranchZPt(DiMuPlottingSystem* dps)
{
  TTree* newtree = new TTree("newtree", "new tree");
  newtree->AddFriend(dps->treename, dps->infilename);

  TBranch *newBranch = newtree->Branch("recoCandMassPrime",&dps->recoCandMassPrime,"recoCandMassPrime/F");

  Int_t nentries = dps->tree->GetEntries();
  std::cout << "nentries: " << nentries << std::endl;
  for (Int_t i = 0; i < nentries; i++)
    {
      dps->tree->GetEntry(i);
      if(dps->recoCandPt > 20 && dps->recoCandPt < 60) dps->recoCandMassPrime = dps->recoCandMass + 0.0000917686*(dps->recoCandPt*dps->recoCandPt)*TMath::Sin(-6.43015*TMath::Sqrt(dps->recoCandPt) + 0.210496);
      else dps->recoCandMassPrime = dps->recoCandMass;
      newBranch->Fill();
    }
  dps->tree = newtree;
}

void addMCWeight(TString filename, TString treename)
{

    using namespace std;

    cout << "adding MCWeight to:"<< filename << endl;

    double scale_w;
    float mcWeight;
    int npv;

    TFile *file = new TFile(filename,"UPDATE");  
    TTree *oldtree = (TTree*)file->Get(treename);

    if(oldtree==NULL)
    {
        cout << "Could not find tree " << treeDir << "/" << treename << endl
             << "in file " << file->GetName() << endl;
        return;
    }
  
    oldtree->SetBranchAddress("scale_w",&scale_w);
    oldtree->SetBranchAddress("mcWeight",&mcWeight);
    oldtree->SetBranchAddress("npv",&npv);

    double scaleMC_w = 1.0;

    TBranch *branch = oldtree->Branch("scaleMC_w",&scaleMC_w,"scaleMC_w/D");

    for(int i = 0; i < oldtree->GetEntries(); i++)
    {
        oldtree->GetEntry(i);
      
        double w_npv = (3.57041 + -1.49846*npv + 0.515829*npv*npv + -0.0839209*npv*npv*npv + 0.00719964*npv*npv*npv*npv + -0.000354548*npv*npv*npv*npv*npv + 1.01544e-05*npv*npv*npv*npv*npv*npv + -1.58019e-07*npv*npv*npv*npv*npv*npv*npv + 1.03663e-09*npv*npv*npv*npv*npv*npv*npv*npv);
        
        scaleMC_w = scale_w*mcWeight*w_npv;
        branch->Fill();
    }

    file->cd();
    oldtree->CloneTree()->Write(treename, TObject::kOverwrite);
    file->Close();
  
}

void addMCWeight_data(TString filename, TString treename)
{

    using namespace std;
    
    cout << "Data adding MCWeight to:"<< filename << endl;

    double scale_w;
    float mcWeight;

    TFile *file = new TFile(filename,"UPDATE");  
    TTree *oldtree = (TTree*)file->Get(treename);

    if(oldtree==NULL)
    {
        cout << "Could not find tree " << treeDir << "/" << treename << endl
             << "in file " << file->GetName() << endl;
        return;
    }
  
    oldtree->SetBranchAddress("scale_w",&scale_w);
    oldtree->SetBranchAddress("mcWeight",&mcWeight);

    double scaleMC_w = 1.0;

    TBranch *branch = oldtree->Branch("scaleMC_w",&scaleMC_w,"scaleMC_w/D");

    for(int i = 0; i < oldtree->GetEntries(); i++)
    {
        oldtree->GetEntry(i);      
        
        scaleMC_w = scale_w*mcWeight;
        branch->Fill();
    }

    file->cd();
    oldtree->CloneTree()->Write(treename, TObject::kOverwrite);
    file->Close();
  
}

void addBranch() {
  //TFile *f = TFile::Open("TTbar_madgraphMLM_1000pb_weighted.root","update");
  //TFile *f = TFile::Open("QCD_1000pb_weighted.root","update");
  TFile *f = TFile::Open("GluGluHToBB_M125_13TeV_powheg_pythia8_1000pb_weighted.root","update"); 
  TTree *T = (TTree*)f->Get("otree");
  double AK8Puppijet0_tau21,AK8Puppijet0_msd,AK8Puppijet0_pt;
  double AK8Puppijet0_tau21DDT = 999;
  TBranch *bpt = T->Branch("AK8Puppijet0_tau21DDT",&AK8Puppijet0_tau21DDT,"AK8Puppijet0_tau21DDT/D");
  T->SetBranchAddress("AK8Puppijet0_tau21",&AK8Puppijet0_tau21);
  T->SetBranchAddress("AK8Puppijet0_msd",&AK8Puppijet0_msd);
  T->SetBranchAddress("AK8Puppijet0_pt",&AK8Puppijet0_pt);
  Long64_t nentries = T->GetEntries();
  for (Long64_t i=0;i<nentries;i++) {
     T->GetEntry(i);
     if (AK8Puppijet0_pt > 0. && AK8Puppijet0_msd > 0.)
       AK8Puppijet0_tau21DDT = AK8Puppijet0_tau21 + 0.063*TMath::Log(AK8Puppijet0_msd*AK8Puppijet0_msd/AK8Puppijet0_pt);
     bpt->Fill();
  }
  T->Print();
  T->Write();
  delete f;
}

void updatePbPbBtriggerweight(TString filename, vector<float> w)
{
  auto f = new TFile(filename,"update");

  auto nt = (TTree *)f->Get("nt");

  float csv60, csv80;
  float triggermatched;
  float weight;
  TBranch *bw;

  bw =  nt->Branch("weight",&weight);

  nt->SetBranchAddress("hltCSV60",&csv60);
  nt->SetBranchAddress("hltCSV80",&csv80);
  nt->SetBranchAddress("triggermatched",&triggermatched);
  
  int n = nt->GetEntries();
  int onep = n/100;
  for (int i=0;i<n;i++) {
    if (i%onep==0) cout<<i/onep<<endl;
    nt->GetEntry(i);


    weight = 0;
    if (triggermatched && csv80) weight = w[1];
    if (triggermatched && csv60 && !csv80) weight = w[0];

    bw->Fill();
  }

  nt->Write();
  f->Close();


}

//.........这里部分代码省略.........
      var8 = intVar8;
      var9 = intVar9;
      var10 = dVar10;
      var11 = dVar11;
      var12 = dVar12;
      var13 = dVar13;
      var14 = dVar14;
      var15 = dVar15;
      if(vecVar16->size()>0){
      	var16 = vecVar16->at(0);
      }
      else{
      	var16 = 0;
      }
      if(vecVar17->size()>0){
      	var17 = vecVar17->at(0);
      }
      else{
      	var18 = 0;
      }
      if(vecVar19->size()>0){
      	var19 = vecVar19->at(0);
      }
      else{
      	var19 = 0;
      }     
      if(vecVar20->size()>0){
      	var20 = vecVar20->at(0);
      }
      else{
      	var20 = 0;
      }
      if(vecVar16->size()>1){
	  	var21 = vecVar16->at(1);
	  }
	  else{
	  	var21 = 0;
	  }
	  var22 = dVar22;
	  var23 = dVar23;
	  var24 = dVar24;
	  var25 = dVar25;
	  var26 = dVar26;
	  var27 = dVar27;
	  var28 = dVar28;
	  var29 = dVar29;      // --- Return the MVA outputs and fill into histograms

      if (Use["CutsGA"]) {
         // Cuts is a special case: give the desired signal efficienciy
         Bool_t passed = reader->EvaluateMVA( "CutsGA method", effS );
         if (passed) nSelCutsGA++;
      }
      BDT = reader->EvaluateMVA( "BDT method");
      histBdt->Fill(BDT);
      branchBDT->Fill();
   }

   // Get elapsed time
   sw.Stop();
   std::cout << "--- End of event loop: "; sw.Print();

   // Get efficiency for cuts classifier
   if (Use["CutsGA"]) std::cout << "--- Efficiency for CutsGA method: " << double(nSelCutsGA)/theTree->GetEntries()
                                << " (for a required signal efficiency of " << effS << ")" << std::endl;

   if (Use["CutsGA"]) {

      // test: retrieve cuts for particular signal efficiency
      // CINT ignores dynamic_casts so we have to use a cuts-secific Reader function to acces the pointer  
      TMVA::MethodCuts* mcuts = reader->FindCutsMVA( "CutsGA method" ) ;

      if (mcuts) {      
         std::vector<Double_t> cutsMin;
         std::vector<Double_t> cutsMax;
         mcuts->GetCuts( 0.7, cutsMin, cutsMax );
         std::cout << "--- -------------------------------------------------------------" << std::endl;
         std::cout << "--- Retrieve cut values for signal efficiency of 0.7 from Reader" << std::endl;
         for (UInt_t ivar=0; ivar<cutsMin.size(); ivar++) {
            std::cout << "... Cut: " 
                      << cutsMin[ivar] 
                      << " < \"" 
                      << mcuts->GetInputVar(ivar)
                      << "\" <= " 
                      << cutsMax[ivar] << std::endl;
         }
         std::cout << "--- -------------------------------------------------------------" << std::endl;
      }
   }

   // --- Write histograms

   newTree->Write("",TObject::kOverwrite);
   target->Close();

   std::cout << "--- Created root file: \""<<oFileName<<"\" containing the MVA output histograms" << std::endl;
  
   delete reader;
    
   std::cout << "==> TMVAClassificationApplication is done!" << endl << std::endl;
} 

//.........这里部分代码省略.........
 int numberOfSamples = ReadFileWeight(nameFileIn, nameSample, PUScenarioIndex); 
 
 std::cout << std::endl;
 std::cout << " numberOfSamples = " << numberOfSamples << std::endl;
 std::cout << std::endl;
 
//  std::cout << " PUScenarioIndex[0] = " << PUScenarioIndex[0]  << std::endl;
   
 ///==== PU reweight (begin) ====
 int numPUScenarios   = gConfigParser -> readIntOption("PU::PUScenarios");

 std::vector<std::string> PUVAR;
 std::vector<double> PUMC[100];
 std::vector<double> PUDATA[100];
 PUclass PU[100];
 
 for (int iScPU = 0; iScPU < numPUScenarios; iScPU++) {
  TString namePU;
  namePU = Form("PU::PUMC%d",iScPU);
  PUMC[iScPU]   = gConfigParser -> readDoubleListOption(namePU.Data());
  namePU = Form("PU::PUDATA%d",iScPU);
  PUDATA[iScPU] = gConfigParser -> readDoubleListOption(namePU.Data());
  namePU = Form("PU::PUVAR%d",iScPU);
  std::string tempName = gConfigParser -> readStringOption(namePU.Data());
  PUVAR.push_back(tempName);
  
  if (PUMC[iScPU].size() != PUDATA[iScPU].size()) {
   std::cout << " PUVAR["  << iScPU << "]        = " << PUVAR.at(iScPU) << std::endl;
   std::cout << " PUMC["   << iScPU << "].size() = " << PUMC[iScPU].size()   << std::endl;
   std::cout << " PUDATA[" << iScPU << "].size() = " << PUDATA[iScPU].size() << std::endl;
   std::cerr << " ERROR " << std::endl;
   return 1;
  }
  
  double sumPUMC = 0;
  for (int itVPU = 0; itVPU < PUMC[iScPU].size(); itVPU++ ){
   sumPUMC += PUMC[iScPU].at(itVPU);  
  }
  double sumPUDATA = 0;
  for (int itVPU = 0; itVPU < PUDATA[iScPU].size(); itVPU++ ){
   sumPUDATA += PUDATA[iScPU].at(itVPU);  
  } 
  
  for (int itVPU = 0; itVPU < PUMC[iScPU].size(); itVPU++ ){
   PU[iScPU].PUWeight.push_back(PUDATA[iScPU].at(itVPU) / PUMC[iScPU].at(itVPU) * sumPUMC / sumPUDATA);
  }
 }
 
 ///==== PU reweight (end) ====
 
 
 ///==== debug flag ==== 
 bool  debug = false; 
 try {
  debug = gConfigParser -> readBoolOption("Input::debug");
 }
 catch (char const* exceptionString){
  std::cerr << " exception = " << exceptionString << std::endl;
 }
 std::cout << ">>>>> input::debug  " << debug  << std::endl;  
 
 
 for (int iSample=0; iSample<numberOfSamples; iSample++){
  std::cout << " nameSample[" << iSample << ":" << numberOfSamples << "] = " << nameSample[iSample] << std::endl;
  
  ///==== GetTree (begin) ==== 
  char nameFile[20000];
  sprintf(nameFile,"%s/out_%s.root",inputDirectory.c_str(),nameSample[iSample]);  
  std::cout << " nameFile = " << nameFile << std::endl;
  TFile* f = new TFile(nameFile, "update");
  f->cd(treeNameDir.c_str());
  
  treeJetLepVect[iSample] = (TTree*) f->Get(treeName.c_str());
//   char nameTreeJetLep[100];
//   sprintf(nameTreeJetLep,"treeJetLep_%d",iSample); 
//   treeJetLepVect[iSample]->SetName(nameTreeJetLep);
  ///==== GetTree (end) ====
  
  std::cout << " >>>>> GetEntries " << treeJetLepVect[iSample]->GetEntries() << std::endl;   
  
  double weight;
//   TBranch *newBranch;
  ///==== add new branch ====
  TBranch *newBranch = treeJetLepVect[iSample] -> Branch("weight_PU",&weight,"weight_PU/D");
  int numPUMC;
  treeJetLepVect[iSample] -> SetBranchAddress(PUVAR.at(PUScenarioIndex[iSample]).c_str(),&numPUMC);
  Long64_t nentries = treeJetLepVect[iSample]->GetEntries();
  for (Long64_t iEntry = 0; iEntry < nentries; iEntry++){
   if((iEntry%((nentries+10)/10)) == 0) std::cout << ">>>>> analysis::GetEntry " << iEntry << " : " << nentries << std::endl;   
   treeJetLepVect[iSample]->GetEntry(iEntry);
   weight = PU[ PUScenarioIndex[iSample] ].getPUWeight(numPUMC);
   newBranch->Fill();
  }
  // save only the new version of the tree
  treeJetLepVect[iSample]->Write("", TObject::kOverwrite);
 }
 
 
 std::cout << " *** end *** " << std::endl;
}

void flagMultCands( TString fname, TString tname ) {

  TFile *inFile = new TFile( fname, "UPDATE" );
  TTree *tree = (TTree*)inFile->Get( tname );

  if ( tree->FindBranch( "pass_multcand" ) ) {
    cout << "pass_multcand branch already exists so I won't add it" << endl;
    delete inFile;
    return;
  }

  ULong64_t  eventNumber;
  int        itype;
  bool       pass_bdt;
  bool       pass_pid;
  bool       pass_rhokst;
  bool       pass_massveto;
  UInt_t     nCandidate;
  ULong64_t  totCandidates;
  tree->SetBranchAddress(  "eventNumber"                 , &eventNumber                 );
  tree->SetBranchAddress(  "itype"                       , &itype                       );
  tree->SetBranchAddress(  "pass_bdt"                    , &pass_bdt                    );
  tree->SetBranchAddress(  "pass_pid"                    , &pass_pid                    );
  tree->SetBranchAddress(  "pass_rhokst"                 , &pass_rhokst                 );
  tree->SetBranchAddress(  "pass_massveto"               , &pass_massveto               );
  tree->SetBranchAddress(  "totCandidates"               , &totCandidates               );
  tree->SetBranchAddress(  "nCandidate"                  , &nCandidate                  );

  map< ULong64_t, ULong64_t > multCandEventNumbers;

  cout << "Finding multiple candidates" << endl;
  // first loop tree and save number of all multiple candidates
  for ( int ev=0; ev<tree->GetEntries(); ev++ ) {
    tree->GetEntry(ev);
    if ( ev%10000==0 ) cout << ev << "/" << tree->GetEntries() << endl;
    //if ( itype>0 && pass_bdt && pass_pid && (!pass_rhokst) && (!pass_massveto) && totCandidates > 1 ) {
    if ( itype>0 && pass_bdt && pass_pid && (!pass_massveto) && totCandidates > 1 ) {
      multCandEventNumbers[eventNumber] = totCandidates;
    }
  }

  // now randomly select which one to keep
  map< ULong64_t, UInt_t > eventToKeep;
  TRandom3 rand;
  rand.SetSeed(2016);
  for ( map<ULong64_t,ULong64_t>::iterator it=multCandEventNumbers.begin(); it!=multCandEventNumbers.end(); it++) {
    UInt_t keep = rand.Integer( it->second );
    eventToKeep[it->first] = keep;
  }

  cout << "Adding flag pass_multcand to tree" << endl;
  // then loop again and write the random choice back in
  bool pass_multcand;
  TBranch *bpt = tree->Branch( "pass_multcand", &pass_multcand, "pass_multcand/O" );

  for ( int ev=0; ev<tree->GetEntries(); ev++ ) {
    tree->GetEntry(ev);
    if ( ev%10000==0 ) cout << ev << "/" << tree->GetEntries() << endl;
    pass_multcand = true;
    if ( totCandidates > 1 && eventToKeep[eventNumber] != nCandidate ) pass_multcand = false;
    bpt->Fill();
  }
  tree->Write();
  delete inFile;
}

//.........这里部分代码省略.........
   Float_t BDT;
   TBranch *branchBDT = newTree->Branch("BDT_"+bkgSample,&BDT,"BDT/F");
   std::vector<Double_t> *vecVar1;
   std::vector<Double_t> *vecVar5;
   std::vector<Double_t> *vecVar7;
   theTree->SetBranchAddress( "svMass", &vecVar1);
   theTree->SetBranchAddress( "dRTauTau", &var3);
   theTree->SetBranchAddress( "dRJJ", &var4 );
//    theTree->SetBranchAddress( "svPt", &vecVar5 );
//    theTree->SetBranchAddress( "dRhh", &var6 );
   theTree->SetBranchAddress( "met", &vecVar7 );
   theTree->SetBranchAddress( "mJJ", &var8 );
//    theTree->SetBranchAddress( "metTau1DPhi", &var9 );
//    theTree->SetBranchAddress( "metTau2DPhi", &var10);
//    theTree->SetBranchAddress( "metJ1DPhi", &var11);
//    theTree->SetBranchAddress( "metJ2DPhi", &var12 );
//    theTree->SetBranchAddress( "metTauPairDPhi", &var13 );
//    theTree->SetBranchAddress( "metSvTauPairDPhi", &var14 );
//    theTree->SetBranchAddress( "metJetPairDPhi", &var15 );
//    theTree->SetBranchAddress( "CSVJ1", &var16 );
//    theTree->SetBranchAddress( "CSVJ2", &var17 );
   theTree->SetBranchAddress( "fMassKinFit", &var2);
   theTree->SetBranchAddress( "chi2KinFit2", &var18);

   //to get initial pre-processed events
   TH1F* cutFlow = (TH1F*)input->Get("preselection");

   // Efficiency calculator for cut method
   Int_t    nSelCutsGA = 0;
   Double_t effS       = 0.7;

   std::vector<Float_t> vecVar(4); // vector for EvaluateMVA tests

   std::cout << "--- Processing: " << theTree->GetEntries() << " events" << std::endl;
   TStopwatch sw;
   sw.Start();
   for (Long64_t ievt=0; ievt<theTree->GetEntries();ievt++) {

      if (ievt%1000 == 0) std::cout << "--- ... Processing event: " << ievt << std::endl;
      theTree->GetEntry(ievt);
      var1 = vecVar1->at(0);
//       var5 = vecVar5->at(0);
      var7 = vecVar7->at(0);
      // --- Return the MVA outputs and fill into histograms

      if (Use["CutsGA"]) {
         // Cuts is a special case: give the desired signal efficienciy
         Bool_t passed = reader->EvaluateMVA( "CutsGA method", effS );
         if (passed) nSelCutsGA++;
      }
      BDT = reader->EvaluateMVA( "BDT method");
      histBdt->Fill(BDT);
      branchBDT->Fill();
   }

   // Get elapsed time
   sw.Stop();
   std::cout << "--- End of event loop: "; sw.Print();

   // Get efficiency for cuts classifier
   if (Use["CutsGA"]) std::cout << "--- Efficiency for CutsGA method: " << double(nSelCutsGA)/theTree->GetEntries()
                                << " (for a required signal efficiency of " << effS << ")" << std::endl;

   if (Use["CutsGA"]) {

      // test: retrieve cuts for particular signal efficiency
      // CINT ignores dynamic_casts so we have to use a cuts-secific Reader function to acces the pointer  
      TMVA::MethodCuts* mcuts = reader->FindCutsMVA( "CutsGA method" ) ;

      if (mcuts) {      
         std::vector<Double_t> cutsMin;
         std::vector<Double_t> cutsMax;
         mcuts->GetCuts( 0.7, cutsMin, cutsMax );
         std::cout << "--- -------------------------------------------------------------" << std::endl;
         std::cout << "--- Retrieve cut values for signal efficiency of 0.7 from Reader" << std::endl;
         for (UInt_t ivar=0; ivar<cutsMin.size(); ivar++) {
            std::cout << "... Cut: " 
                      << cutsMin[ivar] 
                      << " < \"" 
                      << mcuts->GetInputVar(ivar)
                      << "\" <= " 
                      << cutsMax[ivar] << std::endl;
         }
         std::cout << "--- -------------------------------------------------------------" << std::endl;
      }
   }

   // --- Write histograms

   histBdt->Write();
   cutFlow->Write();
   newTree->Write();
   target->Close();

   std::cout << "--- Created root file: \""<<oFileName<<"\" containing the MVA output histograms" << std::endl;
  
   delete reader;
    
   std::cout << "==> TMVAClassificationApplication is done!" << endl << std::endl;
} 

//.........这里部分代码省略.........
    ss >> dummy;
    valHeight.push_back(dummy);
  }

  //Initialise the output tree.
  TChain* OutChain = new TChain("ChainResults", "The collected run results");

  for(size_t i=0; i<theEnergy.size(); i++){
    for(size_t j=0; j<theHeight.size(); j++){

      string fname = "Edep_" + theEnergy.at(i) + "MeV_" + theHeight.at(j)
	+ "cm.root";

      cout << "Processing File: " << fname << "... ";

      TFile* infile = TFile::Open(fname.c_str(), "update");
      TTree* inTree = (TTree*)infile->Get("Results");

      //Clean up any old branches...
      //TBranch* OldBranch;
      //if(inTree->GetListOfBranches()->FindObject("BeamEn")){
      //OldBranch = inTree->GetBranch("BeamEn");
      //inTree->GetListOfBranches()->Remove(OldBranch);
      //inTree->Write();
      //}
      //if(inTree->GetListOfBranches()->FindObject("BeamHgt")){
      //OldBranch = inTree->GetBranch("BeamHgt");
      //inTree->GetListOfBranches()->Remove(OldBranch);
      //inTree->Write();
      //}
      //if(inTree->GetListOfBranches()->FindObject("MeanPE_T1")){
      //OldBranch = inTree->GetBranch("MeanPE_T1");
      //inTree->GetListOfBranches()->Remove(OldBranch);
      //inTree->Write();
      //}
      //if(inTree->GetListOfBranches()->FindObject("MeanPE_T2")){
      //OldBranch = inTree->GetBranch("MeanPE_T2");
      //inTree->GetListOfBranches()->Remove(OldBranch);
      //inTree->Write();
      //}
      

      float BeamEn = 0;
      float BeamHgt = 0;
      float MeanT1 = 0;
      float MeanT2 = 0;
      //Add an energy and beam height branch to the tree.
      TBranch* EnBranch = inTree->Branch("BeamEn", &BeamEn, "BeamEn/F");
      TBranch* HgtBranch = inTree->Branch("BeamHgt", &BeamHgt, "BeamHgt/F");
      TBranch* MeanPE_T1 = 
	inTree->Branch("MeanPE_T1", &MeanT1, "MeanPE_T1/F");
      TBranch* MeanPE_T2 =
	inTree->Branch("MeanPE_T2", &MeanT2, "MeanPE_T2/F");
      unsigned long long int measPE_T1 = 0;
      unsigned long long int measPE_T2 = 0;
      double Edep_H1 =0;
      double Edep_H2 = 0;

      inTree->SetBranchAddress("MeasNumPhotons_PMTT1", &measPE_T1);
      inTree->SetBranchAddress("MeasNumPhotons_PMTT2", &measPE_T2);
      inTree->SetBranchAddress("Edep_H1", &Edep_H1);      
      inTree->SetBranchAddress("Edep_H2", &Edep_H2);
      Long64_t nentries = inTree->GetEntries();
      int numAvg = 0;
      for(int k = 0; k<nentries; k++){
	inTree->GetEntry(k);
	//Calculate means
	if((Edep_H1!=0)&&(Edep_H2!=0)){
	  MeanT1 += measPE_T1;
	  MeanT2 += measPE_T2;
	  numAvg++;
	}
      }

      MeanT1 = MeanT1/numAvg;
      MeanT2 = MeanT2/numAvg;

      for(int k = 0; k<nentries; k++){
	BeamEn = valEnergy.at(i);
	BeamHgt = valHeight.at(j);
	EnBranch->Fill();
	HgtBranch->Fill();
	MeanPE_T1->Fill();
	MeanPE_T2->Fill();
      }

      inTree->Write("", TObject::kOverwrite);
	
      fname = fname + "/Results";
      OutChain->Add(fname.c_str());
      infile->Close();
      cout << "Done!" << endl;
    }
  }
  cout << "OutChain has " << OutChain->GetEntries() << " entries." << endl;
  TFile* OutFile = new TFile("ChainedRuns.root", "RECREATE");
  OutChain->Write();
  OutFile->Close();
  return;
}

//.........这里部分代码省略.........
            //const AliPHOSDigit *digit = prl->Digit(iDig);
            AliPHOSDigit *digit = (AliPHOSDigit*)phosDigits->At(iDig);
            nDigits++;
            for (Int_t n = 0; n < nMyDigits; n++)
            {
                AliPHOSDigit *myDigit = (AliPHOSDigit*)mydigits->At(n);
                if (digit->GetId() == myDigit->GetId())
                {
                    nOverlappingDigits++;
                    break;
                }
            }
        }
        if(nOverlappingDigits == nMyDigits)
        {
            std::cout << "Digits alredy embedded!" << std::endl;
            continue;
        }
        for (Int_t iDig = 0; iDig < nMyDigits; iDig++)
        {
            AliPHOSDigit *myDigit = (AliPHOSDigit*)mydigits->At(iDig);
            if (myDigit)
            {
                for (Int_t n = 0; n < nPhosDigits; n++)
                {
                    //const AliPHOSDigit *digit = prl->Digit(n);

                    AliPHOSDigit *digit = (AliPHOSDigit*)phosDigits->At(n);
                    if (digit->GetId() == myDigit->GetId())
                    {
                        digit->SetALTROSamplesHG(0, 0);
                        digit->SetALTROSamplesLG(0, 0);

                        *digit += *myDigit;
                        myDigit = 0;
                        break;
                    }
                }
                if (myDigit)
                {
                    TClonesArray *digArray = prl->Digits();
                    AliPHOSDigit *newDig =  new((*digArray)[nPhosDigits+nNewDigits]) AliPHOSDigit(*myDigit);

                    newDig->SetALTROSamplesHG(0, 0);
                    newDig->SetALTROSamplesLG(0, 0);
                    nNewDigits++;
                }
            }
        }
        phosDigits->Compress();
        Int_t ndigits = phosDigits->GetEntries() ;
	phosDigits->Sort();
        // Remove digits that are flagged bad in BCM. Then remove digits that are below threshold
        for (Int_t i = 0 ; i < ndigits ; i++)
        {
            AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( phosDigits->At(i) ) ;
	    //std::cout << digit->GetId() << std::endl;
            if(digit->GetId())
            {
                vector<Int_t>::iterator it;
                it = std::find (badChannels.begin(), badChannels.end(), digit->GetId() );
                if(*it)
                {
                    digit->SetEnergy(0.0);
                }
            }
            if(digit->GetEnergy() <= recoParam->GetGlobalAltroThreshold())
	    {
	      phosDigits->RemoveAt(i);
	    }
        }
        //Set indexes in list of digits and make true digitization of the energy
        phosDigits->Compress();
        phosDigits->Sort();
	ndigits = phosDigits->GetEntries();
        for (Int_t i = 0 ; i < ndigits ; i++)
        {
            AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( phosDigits->At(i) ) ;
            digit->SetIndexInList(i) ;
        }
        // -- create Digits branch
        Int_t bufferSize = 32000 ;

        TObjArray *branchList = prl->TreeD()->GetListOfBranches();

        branchList->RemoveAt(0);

        TBranch * digitsBranch = prl->TreeD()->Branch("PHOS","TClonesArray",&phosDigits,bufferSize);

        digitsBranch->Fill() ;
        prl->WriteDigits("OVERWRITE");
    }
    prl->WriteDigits("OVERWRITE");
    std::cout << "# Digits: " << nDigits << std::endl;
    std::cout << "# Embedded digits: " << nEmbedDigits << std::endl;
    std::cout << "# Overlapping digits: " << nOverlappingDigits << std::endl;
    std::cout << "# New digits: " << nNewDigits << std::endl;
    
    return 0;
}