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

void fakeDAQ(const char* outFile = "/Users/chaozhang/Projects/LBNE/WbLS/software/WbLSDAQ/data/test/fakedata.root")
{
    // WblsDaq::Spinner spinner("/Users/chaozhang/Projects/LBNE/WbLS/software/WbLSDAQ/data/test/example2.root");
    // TTree* header = spinner.header_tree();
    // TTree* footer = spinner.footer_tree();
    // TTree* events = spinner.event_tree();
    
    TFile* fp = TFile::Open("/Users/chaozhang/Projects/LBNE/WbLS/software/WbLSDAQ/data/test/example2.root"); // intentional leak
    TTree* header = (TTree*)fp->Get("Header");
    TTree* footer = (TTree*)fp->Get("Footer");
    TTree* events = (TTree*)fp->Get("FADCData");
    
    const int N = 120; // fake 400 events
    int nentries = events->GetEntries();
    gRandom->SetSeed(0);
    int startN = int(gRandom->Uniform(nentries));
    TString sel = Form("Entry$>=%d && Entry$<%d", startN, startN + N);
    TFile f(outFile, "recreate");    
    TTree *th = header->CopyTree("Entry$<8");
    TTree *tf = footer->CopyTree("Entry$<8");
    TTree *te = events->CopyTree(sel.Data());
    // TTree *te = events->CopyTree("Entry$>1850 && Entry$<1970");
    // cout << sel << endl;
    th->Write();
    tf->Write();
    te->Write();
    f.Close();
    // cout << outFile << endl;
}

void Skimmer(TString cname)
{
  TFile* file=TFile::Open(cname + "/treeProducerStop4Body/tree.root", "read");
  TString cutString = "nLepGood >= 1"; // Place custom cut string for the skim here

  if(!file)
  {
    std::cout << "Corrupted file! : " << cname << std::endl;
    return;
  }

  TTree* tree = static_cast<TTree*>(file->Get("tree"));
  if(!tree)
  {
    std::cout << "Corrupted file! : " << cname << std::endl;
    return;
  }

  TFile* ofile = new TFile(cname + "/treeProducerStop4Body/tree_skimmed.root", "recreate");
  TTree* otree = tree->CopyTree(cutString);

  ofile->Write();
  ofile->Close();

  gROOT->ProcessLine(".q");
}

void skimFastNtuple(TString inputFile = "input.root",TString outputFile = "output.root", int skim = 0){

TFile *_file0 = TFile::Open(inputFile);
_file0->cd();
TTree* tree = (TTree*)_file0->Get("tree");
TTree* oldtree = tree;
cout << oldtree->GetEntries() << endl;

TFile *_file1 = TFile::Open(outputFile,"recreate");
_file1->cd();

char cutString[200];
if     (skim == 0) sprintf(cutString,"met>20||lid3!=0");
else if(skim == 1) sprintf(cutString,"(type==1||type==2)&&lid3==0&&njets<=1");
else if(skim == 2) sprintf(cutString,"min(pmet,pTrackMet)>20&&(TMath::Abs(dilep->mass()-91.1876)>15||type==1||type==2)&&lid3==0&&dilep->pt()>45&&njets<=3");
else if(skim == 3) sprintf(cutString,"((cuts & 4) == 4) && ((cuts & 512) == 512)");
else if(skim == 4) sprintf(cutString,"((cuts & 4) != 4) || ((cuts & 512) != 512)");
else if(skim == 5) sprintf(cutString,"lid3!=0");
else if(skim == 6) sprintf(cutString,"min(pmet,pTrackMet)>20&&((TMath::Abs(dilep->mass()-91.1876)>15&&dilep->pt()>45)||type==1||type==2)&&lid3==0&&dilep->pt()>30&&njets<=3");
else if(skim == 7) sprintf(cutString,"((cuts & 4) != 4) && ((cuts & 512) != 512)");
else if(skim == 8) sprintf(cutString,"lq1*lq2>0&&lid3==0");
else if(skim == 9) sprintf(cutString,"min(pmet,pTrackMet)>20&&((TMath::Abs(dilep->mass()-91.1876)>15&&((dymva>0.88&&njets==0)||(dymva>0.84&&njets==1)||(met>45&&njets>=2)))||type==1||type==2)&&lid3==0&&dilep->M()>12&&dilep->pt()>0");
else if(skim ==10) sprintf(cutString,"dilep->pt()>70&&TMath::Abs(dilep->mass()-91.1876)<30&&njets<=1");
else if(skim ==11) sprintf(cutString,"dilep->pt()>70&&njets<=1");
else {cout << "No good option" << endl; return;}
cout << "cut: " << cutString << endl;

TTree *newtree1;
newtree1 = oldtree->CopyTree(cutString);
cout << newtree1->GetEntries() << endl;

_file1->Write();
_file1->Close();
}

void Copy(){

  TFile* tfin = new TFile("/Volume0/ExpData/2012_Feb_Beam/RootFile_cosmic/CosmicResult_20120209.root");

  TFile* tfOut = new TFile("test.root","Recreate");
  TTree* tr   = (TTree*)tfin->Get("GainFitPar");
  TTree* trClone = tr->CopyTree("");
  trClone->Write();
  tfOut->Close();

}

TTree* read_tree(char* filename, char* tree_name, char* tree_cuts)
{
    DEBUG("now reading tree " << tree_name << " with cuts: " << tree_cuts << " from file: " << filename);
    TFile* file = new TFile(filename,"READ");
    TTree* tree = (TTree*)file -> Get(tree_name);
    TFile* fileTmp = new TFile("tmpFile_.root","RECREATE");
    TTree* smallerTree= tree->CopyTree(tree_cuts);
    file -> Close();
    delete file;
    return smallerTree;
}

void CreateFriendTree() {
   // Open the file created by CreateParentTree
   // Copy a subset of the TTree into a new TTree
   //   (see also tutorials copytree.C, copytree2.C and copytree3.C)
   // Create an index on the new TTree ("Run","Event")
   // Write the new TTree (including its index)
   
   TFile *f = new TFile("treeparent.root");
   TTree *T = (TTree*)f->Get("T");
   TFile *ff = new TFile("treefriend.root","recreate");
   TTree *TF = T->CopyTree("z<10");
   TF->SetName("TF");
   TF->BuildIndex("Run","Event");
   TF->Write();
   TF->Print();
   delete ff;
}

void skimSimple2020_20m50(char* fname){
  gSystem->Load("libFWCoreFWLite");
  AutoLibraryLoader::enable();
  
  TFile* tf = new TFile(fname);
  TTree* tr = (TTree*)tf->Get("Events");

  TFile* tfCp = new TFile(Form("%s_skimSimple2020_20m50",fname), "recreate");

  std::string selection = "Sum$(hyp_ll_p4.pt()>20&&hyp_lt_p4.pt()>20)>0&&Sum$(genps_id==23&&genps_p4.M()>20.&&genps_p4.M()<50.)>0";
  int nEvents = tr->GetEntries(); 
  int skimEvents = tr->GetEntries(selection.c_str()); 
  std::cout<<selection.c_str()<<" will select "<<skimEvents<<" out of "<<nEvents<<std::endl;   

  trCp = tr->CopyTree(selection.c_str(), "");

  trCp->Write();

  tfCp->Write();

  std::cout<<"Closing file "<<tfCp->GetName()<<" with "<<trCp->GetEntries()<<" entries in "<<trCp->GetName()<<std::endl;
  tfCp->Close();
}

void ApplyMoreCutsBeforeFit(){
  //------------------------------------------------
    TFile * f = new TFile ("/exp/LHCb/amhis/LeptonU-/tuples/data/forfit/real-data-electrons-photon-30032016-reduced.root");
  //TFile * f = new TFile ("/exp/LHCb/amhis/LeptonU-/tuples/data/forfit/real-data-electrons-jpsi-30032016-reduced.root");
  TTree * signaltree = (TTree*)f->Get("DecayTree");
  TCut PID = "Kaon_MC12TuneV3_ProbNNk > 0.1 && Proton_MC12TuneV3_ProbNNp>0.1&&L2_PIDe>1.2&&L1_PIDe>1.2";
  
  
  // Cuts taken from  pk mumu analysis
  //  TCut HarderPID= "Proton_MC12TuneV3_ProbNNp >0.2 &&Proton_MC12TuneV3_ProbNNk<0.8 && Proton_MC12TuneV3_ProbNNpi<0.7 && Kaon_MC12TuneV3_ProbNNk>0.2 && Kaon_MC12TuneV3_ProbNNp < 0.8 &&L2_PIDe>1.2&&L1_PIDe>1.2" ;
  

  TCut HarderPID= "Proton_ProbNNp >0.2 &&Proton_ProbNNk<0.8 && Proton_ProbNNpi<0.7 && Kaon_ProbNNk>0.2 && Kaon_ProbNNp < 0.8 &&L2_PIDe>1.2&&L1_PIDe>1.2" ;
  
  TCut Kinematics = "Proton_PT >500 && Lambdab_DTF_PV_chi2<10  && Lambdab_PT>2000 ";

  TCut Lstar = "Lambdastar_M<1550 && Lambdastar_M>1450";
 
  //  TCut HopCut = "Lambdab_MM> (4000 + 344 +  100*log(Lambdab_FDCHI2_OWNPV)) ";

  TCut HopCut = "Lambdab_MM> (4000 + 344 +  75*log(Lambdab_FDCHI2_OWNPV)) ";

  TCut RemovePhi = "Lambdab_M01_Subst0_p2K>1050";
  //  TCut RemoveKst = "Lambdab_M01_Subst0_p2pi< 850 && Lambdab_M01_Subst0_p2pi>920";
  TCut RemoveKst = "Lambdab_M01_Subst0_p2pi< 850 ";
  
  //  TCut ApplyTheseCuts = HarderPID && Kinematics&& HopCut && RemovePhi && RemoveKst;
  TCut ApplyTheseCuts = HarderPID && Kinematics ;
  

  TFile * localFileForFit = new TFile("localFileForFit.root", "RECREATE"); 
  cout <<" These are the cuts that we applied : " <<  ApplyTheseCuts << endl;
  TTree * smallerTreeForFit = signaltree->CopyTree(ApplyTheseCuts);
  smallerTreeForFit->Write();
  f->Close();
  return;
}

void plot(bool sgn_ = true, bool doROC_ = false, int steps_ = 50){


  TCanvas *c1 = new TCanvas("c1","",5,30,650,600);
  c1->SetGrid(0,0);
  c1->SetFillStyle(4000);
  c1->SetFillColor(10);
  c1->SetTicky();
  c1->SetObjectStat(0);

  TLegend* leg = new TLegend(0.55,0.47,0.85,0.85,NULL,"brNDC");
  leg->SetFillStyle(0);
  leg->SetBorderSize(0);
  leg->SetFillColor(10);
  leg->SetTextSize(0.04);

  TFile *fSgn = new TFile("/data_CMS/cms/lbianchini/VbfJetsStudy/OpenNtuples/MuTauStream2011_iter2/Inclusive/nTupleDYJets-Mu-50-madgraph-PUS1_Open_MuTauStream.root","READ");
  TFile *fBkg = new TFile("/data_CMS/cms/lbianchini/VbfJetsStudy/OpenNtuples/MuTauStream2011_iter2/Inclusive/nTupleQCDmu_Open_MuTauStream.root","READ");

  if((fSgn->IsZombie()) || (fBkg->IsZombie()) ){
    cout << "No file" << endl;
    exit(1);
  }

  float vx_combRelIsoLeg1Raw[steps_];
  float vy_combRelIsoLeg1Raw[steps_];
  float vx_combRelIsoLeg1[steps_];
  float vy_combRelIsoLeg1[steps_];
  float vx_combRelIsoLeg1Beta[steps_];
  float vy_combRelIsoLeg1Beta[steps_];
  float vx_combRelIsoLeg1DBeta[steps_];
  float vy_combRelIsoLeg1DBeta[steps_];
  float vx_combRelIsoLeg1Rho[steps_];
  float vy_combRelIsoLeg1Rho[steps_];

  TH2F* hMaster = new TH2F("hMaster"," ; signal efficiency ; background efficiency", 100,0,1,100,0,1);

  TTree* treeSgn = (TTree*)fSgn->Get("outTreePtOrd");
  TTree* treeBkg = (TTree*)fBkg->Get("outTreePtOrd");

  TFile* dummy = new TFile("dummy.root","RECREATE");
  //TTree* treeSgnCut = (TTree*)treeSgn->CopyTree("isTauLegMatched && pt1>25 && pt2>20 && eta1*eta2<0");  
  //TTree* treeBkgCut = (TTree*)treeBkg->CopyTree("tightestHPSWP>-1 && pt1>25 && pt2>20 && eta1*eta2<0");  
  TTree* treeSgnCut = (TTree*)treeSgn->CopyTree("isTauLegMatched");  
  TTree* treeBkgCut = (TTree*)treeBkg->CopyTree("tightestHPSWP>0");  


  cout << treeSgnCut->GetEntries() << endl;
  cout << treeBkgCut->GetEntries() << endl;

  TH1F* SgnEffVsNvtx = new TH1F("SgnEffVsNvtx","; # reco-vtx; efficiency for CombRelIso<0.1",16,1,17);
  TH1F* BkgEffVsNvtx = new TH1F("BkgEffVsNvtx","; # reco-vtx; efficiency for CombRelIso<0.1",16,1,17);

  std::vector<string> isoDep;
  isoDep.push_back("combRelIsoLeg1");
  isoDep.push_back("combRelIsoLeg1Beta");
  isoDep.push_back("combRelIsoLeg1DBeta");
  isoDep.push_back("combRelIsoLeg1Rho");
  isoDep.push_back("combRelIsoLeg1Raw");

  std::map<std::string, TH1F*> aMapSgn;
  aMapSgn["combRelIsoLeg1Raw"] = (TH1F*)SgnEffVsNvtx->Clone("combRelIsoLeg1Raw");
  aMapSgn["combRelIsoLeg1"] = (TH1F*)SgnEffVsNvtx->Clone("combRelIsoLeg1");
  aMapSgn["combRelIsoLeg1Beta"] = (TH1F*)SgnEffVsNvtx->Clone("combRelIsoLeg1Beta");
  aMapSgn["combRelIsoLeg1DBeta"] = (TH1F*)SgnEffVsNvtx->Clone("combRelIsoLeg1DBeta");
  aMapSgn["combRelIsoLeg1Rho"] = (TH1F*)SgnEffVsNvtx->Clone("combRelIsoLeg1Rho");
  std::map<std::string, TH1F*> aMapBkg;
  aMapBkg["combRelIsoLeg1Raw"] = (TH1F*)BkgEffVsNvtx->Clone("combRelIsoLeg1Raw");
  aMapBkg["combRelIsoLeg1"] = (TH1F*)BkgEffVsNvtx->Clone("combRelIsoLeg1");
  aMapBkg["combRelIsoLeg1Beta"] = (TH1F*)BkgEffVsNvtx->Clone("combRelIsoLeg1Beta");
  aMapBkg["combRelIsoLeg1DBeta"] = (TH1F*)BkgEffVsNvtx->Clone("combRelIsoLeg1DBeta");
  aMapBkg["combRelIsoLeg1Rho"] = (TH1F*)BkgEffVsNvtx->Clone("combRelIsoLeg1Rho");

  TH1F* h1 = new TH1F("h1","",1,-10,10);

  if(doROC_){
    for(unsigned int m = 0; m < isoDep.size() ; m++){
    
      treeSgnCut->Draw( "eta1>>h1" ,"puWeight");
      float sgnDen = (float)h1->Integral();
      cout << sgnDen << endl;
      h1->Reset();
      treeBkgCut->Draw( "eta1>>h1" ,"puWeight");
      float bkgDen = (float)h1->Integral();
      cout << bkgDen << endl;
      h1->Reset();
      
      for(int i = 0; i < steps_; i++){
	
      cout << isoDep[m] << " --> " << i << endl;
      
      TCut cut(Form("puWeight*(%s<=%f)", isoDep[m].c_str(), i*(1./steps_)) );
      treeSgnCut->Draw( "eta1>>h1" ,cut );
      float sgnNum = (float)h1->Integral();
      cout << sgnNum << endl;
      h1->Reset();
      treeBkgCut->Draw( "eta1>>h1" ,cut );
      float bkgNum = (float)h1->Integral();
      cout << bkgNum << endl;
      h1->Reset();
//.........这里部分代码省略.........

/*****************************************************************
* root macro to skim MT2trees                                    *
* -> sample and path to shlib given as arguments                 *
*                                                                *
* Pascal Nef                             September 18th, 2011    *  
*****************************************************************/
void MT2treeSkimming(string sample, string shlib, string prefix) {
  	gSystem->Load("libPhysics");
  	gSystem->Load(shlib.c_str());

	string LABEL  = "";
	string file   = sample;
	string outfile = prefix+"/"+sample;

	// log file 
	TString log =sample+".skim.log";
	TString cuts="cuts.skim.log";
	string line="";
	ofstream f_log;
	ofstream f_cuts;
	f_log .open(log.Data(),ios::app);
	f_cuts.open(cuts.Data());
	if(!f_log.is_open()||!f_cuts.is_open()) {cout << "ERROR: cannot open file. " << endl; exit(-1);}

	// cuts --------------------------------------------
  	std::ostringstream cutStream;
	cutStream       << " " 	  
//	  << "misc.MT2 >=50"                                     << "&&" 
//	  << "misc.MET>=30"                                      << "&&"
//	  << "misc.HT > 750 "                                    << "&&"
//	  << "misc.Jet0Pass ==1"                                 << "&&"
//	  << "misc.Jet1Pass ==1"                                 << "&&"
//	  << "misc.SecondJPt  >100"                              << "&&"
//	  << "misc.PassJetID ==1"                                << "&&"
//	  << "misc.Vectorsumpt < 70"                             << "&&"
//	  << "((misc.MinMetJetDPhi >0.3&&NBJets==0)||NBJets>=1)" << "&&"
// Lepton Veto
//	  << "(NEles==0 || ele[0].lv.Pt()<10)"                   << "&&"
//	  << "(NMuons==0 || muo[0].lv.Pt()<10)"                  << "&&"
// Lepton Skim
//	  << "(ele[0].lv.Pt()>10 || muo[0].lv.Pt()>10)"          << "&&"
// LowMT2 ----------------------------
//	  << "misc.LeadingJPt >150"                              << "&&"
//	  << "NBJets >0"                                         << "&&"
//	  << "NJetsIDLoose >=4"                                  << "&&"
// -----------------------------------
//	  << "NJetsIDLoose40 >=3"                                << "&&"
// Photons
	  << "(GenDiLeptPt(0,10,0,1000,false)>=100||GenPhoton[0].Pt()>=100)"  << "&&"
//	  << "(RecoOSDiLeptPt(10,10,0,10000)>=100 ||photon[0].lv.Pt()>=100)"  << "&&"
//	  << "NPhotons >0"  
// Noise
//	  << "misc.HBHENoiseFlag == 0"                           << "&&"
//	  << "misc.CSCTightHaloID==0"                            << "&&"
	  << "misc.CrazyHCAL==0";

	
	TString basecut = cutStream.str();
	string  SEL= "("+basecut+")";
	cout << "Skimming with sel: " << SEL << endl;
	TString cuts_log = basecut.ReplaceAll("&&", "\n");
	f_cuts << "Cuts applied: " <<  cuts_log << endl;
	//  --------------------------------------------


	// files ---------------------------------------
	f_log << "Skimming file: " << file;
	TFile *_file0 = TFile::Open( (file).c_str()); 
	TTree * t = (TTree*) _file0->Get("MassTree");
	TH1F*   h_PUWeights = (TH1F*) _file0->Get("h_PUWeights");
	TH1F*   h_Events    = (TH1F*) _file0->Get("h_Events");
	t->SetMaxTreeSize(19000000000);
	TFile*out = TFile::Open( (outfile).c_str(),"RECREATE");
	TTree *tc = t->CopyTree(SEL.c_str());   
	int nentries = tc->GetEntries();
	f_log << " -> skimmed tree has " << nentries << " entries." <<endl;
	f_log.close();
	f_cuts.close();
	out->Write();
	if(h_PUWeights!=0){
		cout << "writing TH1F h_PUWeights" << endl;
		h_PUWeights->Write();
	}
	if(h_Events!=0){
		cout << "writing TH1F h_Events" << endl;
		h_Events->Write();
	}
	out->Close();
	_file0->Close();
	cout << "Result file: " << outfile << endl;
	// -------------------------------------------------

}

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

  for (Long64_t i=0; i<nentries;i++) {
    if (jobid==0 && i % oneperc == 0 && i>0) {
      std::cout << std::fixed;
      TTimeStamp t1; cout<<" \r"<<i/oneperc<<"%   "<<" est. time "<<setprecision(2) <<(t1-t0)*nentries/(i+.1)<<" s ";
      cout<<";Processing:"<<setprecision(2)<<processingTime/(t1-t0)*100<<" %";
      cout<<";Copy1:"<<setprecision(2) <<copyToTime/(t1-t0)*100<<" %";
      cout<<";Clone:"<<setprecision(2) <<cloneTime/(t1-t0)*100<<" %";
      cout<<";Copy2:"<<setprecision(2) <<copyFromTime/(t1-t0)*100<<" %";
      cout<<";Fill:"<<setprecision(2) <<fillTime/(t1-t0)*100<<" %";
      cout<<";Read:"<<setprecision(2) <<readTime/(t1-t0)*100<<" %";
      cout<<flush; 
    }
    
    TTimeStamp tsRead0;
    tjet->GetEntry(i+nentries1);
    TTimeStamp tsRead1;
    readTime+=tsRead1-tsRead0;

    Everything ev;
    TTimeStamp tsCpTo0;
    for (unsigned j=0;j<brInt.size();j++) ev.PutInt(brInt[j],valIntIn[j]);
    for (unsigned j=0;j<brFloat.size();j++) ev.PutFloat(brFloat[j],valFloatIn[j]);
    for (unsigned j=0;j<brVFloat.size();j++) ev.PutVFloat(brVFloat[j],brVFloatCounter[j],valVFloatIn[j]);
    for (unsigned j=0;j<brVInt.size();j++) ev.PutVInt(brVInt[j],brVIntCounter[j],valVIntIn[j]);
    TTimeStamp tsCpTo1;
    copyToTime+=tsCpTo1-tsCpTo0;
    

    TTimeStamp tsCl0;
    //think about: copy object (timing 10% ->3%) 
    //but it copies vectors, so push_back will add in the end...
    //    Everything evout = ev;
    //or even reference(in place?) (->0.2%)
    //Everything &evout = ev;
    Everything evout = ev.CloneStructure();
    TTimeStamp tsCl1;
    cloneTime+=tsCl1-tsCl0;

    bool exitEvent = false;
    int counter = 0;

    while (!exitEvent) {
    
    TTimeStamp tsPr0;
    if (useOneToOne) {
      fProcessOneToOne(ev, evout);
      evout.UpdateCounters();
      exitEvent = true;
    } else {
      exitEvent = fProcessOneToMany(ev, evout, counter); 
      //      if (!exitEvent) cout<<"event to write "<<counter<<endl;
      counter++;
    }

    TTimeStamp tsPr1;  
    processingTime+=tsPr1-tsPr0;
    
    
    //Everything evout = ev;
    TTimeStamp tsCpFrom0;
    for (unsigned j=0;j<brInt.size();j++) valIntOut[j] = evout.GetInt(brInt[j]);
    for (unsigned j=0;j<brFloat.size();j++) {valFloatOut[j] = evout.GetFloat(brFloat[j]); 
      //  cout<<brFloat[j]<<" "<<evout.GetFloat(brFloat[j])<<endl;
    }
    for (unsigned j=0;j<brVFloat.size();j++) valVFloatOut[j] = evout[brVFloat[j]];
    for (unsigned j=0;j<brVInt.size();j++) valVIntOut[j] = evout.GetVInt(brVInt[j]);
    TTimeStamp tsCpFrom1;
    copyFromTime+=tsCpFrom1-tsCpFrom0;

    TTimeStamp tsFill0;
    tjetout->Fill();
    TTimeStamp tsFill1;
    fillTime+=tsFill1-tsFill0;
  }

  }
  cout<<endl;
  s.Stop();
  cout<<"Done in ";s.Print();
  
  tjetout->FlushBaskets();
  tjetout->Write();
  
  cout<<"Copying other trees..."<<endl;

  for (unsigned i=0;i<friendTrees.size();i++) {
    TTree *t = friendTrees[i];
    if (sameFileFriend[i]) {
      //TTree *triendout = t->CloneTree(-1,"fast");
      TTree *triendout = t->CopyTree("","",nentries,nentries1);
      triendout->Write();
    }
  }
  


  fout->Close();
  friendTrees.clear();
}

int x3d( TString cut, int showFlag )
{

  // Retrieve trees and apply cut
  TFile* alignedFile = new TFile("aligned.root");
  TTree* tmpTree     = (TTree*)alignedFile->Get("theTree");
  TTree* alignedTree = (TTree*)tmpTree->CopyTree(cut);

  TFile* misalignedFile = new TFile("misaligned.root");
  tmpTree        = (TTree*)misalignedFile->Get("theTree");
  TTree* misalignedTree = (TTree*)tmpTree->CopyTree(cut);

  // Set tree branches
  float x,y,z,phi,theta,length,thick,width;
  float mx,my,mz,mphi,mtheta,mlength,mthick,mwidth;
  TRotMatrix* rot;
  TRotMatrix* mrot;
  double rad2deg = 180./3.1415926;

  alignedTree->SetBranchAddress( "x",      &x      );
  alignedTree->SetBranchAddress( "y",      &y      );
  alignedTree->SetBranchAddress( "z",      &z      );
  alignedTree->SetBranchAddress( "phi",    &phi    );
  alignedTree->SetBranchAddress( "theta",  &theta  );
  alignedTree->SetBranchAddress( "length", &length );
  alignedTree->SetBranchAddress( "width",  &width  );
  alignedTree->SetBranchAddress( "thick",  &thick  );
  alignedTree->SetBranchAddress( "rot",    &rot    );

  misalignedTree->SetBranchAddress( "x",      &mx      );
  misalignedTree->SetBranchAddress( "y",      &my      );
  misalignedTree->SetBranchAddress( "z",      &mz      );
  misalignedTree->SetBranchAddress( "phi",    &mphi    );
  misalignedTree->SetBranchAddress( "theta",  &mtheta  );
  misalignedTree->SetBranchAddress( "length", &mlength );
  misalignedTree->SetBranchAddress( "width",  &mwidth  );
  misalignedTree->SetBranchAddress( "thick",  &mthick  );
  misalignedTree->SetBranchAddress( "rot",    &mrot    );

  // Create canvas
  TCanvas* c1 = new TCanvas("c1","Detector units", 200, 10, 700, 500);
  c1->cd();

  TBRIK* IP = new TBRIK("IP","IP","void",0.,0.,0.);
  TNode* rootNode = new TNode("Root","Root","IP",0.,0.,0.);
  rootNode->cd();

  int entry = 0;
  while ( alignedTree->GetEntry(entry) && misalignedTree->GetEntry(entry) )
	{
	  entry++;
	  std::ostringstream name;

	  // Aligned detector
	  name << "aBrik" << entry;
	  TBRIK* aBrik = new TBRIK(name.str().c_str(),"Aligned detector unit","void",
							   0.01,0.01,length);
	  aBrik->SetLineColor(4);

	  // Detector node (position and orientation)
	  name.str("aNode"); name << entry;
	  TNode* aNode = new TNode(name.str().c_str(),name.str().c_str(),aBrik,x,y,z);
	  // Misaligned detector
 	  name.str("mBrik");
 	  name << entry;
 	  TBRIK* mBrik = new TBRIK(name.str().c_str(),"Misaligned detector unit","void",
							   0.01,0.01,mlength);
 	  mBrik->SetLineColor(2);

	  // Detector node (position and orientation)
 	  name.str("mNode"); name << entry;
 	  TNode* mNode = new TNode(name.str().c_str(),name.str().c_str(),mBrik,mx,my,mz);

	  //if (entry>5) break;
	}

  rootNode->cd();
  rootNode->Draw();
  
  c1->GetViewer3D();

  return 0;

}

TTree * InitTrees(const char * detector,  const char *referenceDet){
  //
  // Init tree for given period
  //  all  trees stored in qaMap
  //  FriendTree added to the master treeQADet
  // Currentrly in the code we assume ID="run";
  // Bug in root tree ?  - in case more than one friend tree set - indeces looks corrupted 
  //    0.) QA tree per detector Raw+MC (if exist)
  //    1.) QA trees per refernce detectors specified by string refDet e.g "TPC;TRD;TOF"
  //    2.) Logbook tree per run
  //    3.) Logbook tree per run/detector
  //    3.) RCT table
  //    4.) CPass table 
  // 
  //  tree is created with addition of friend trees which can be used in queries
  //  queries as analog to the SQL statement

  /* 
    period="LHC15o"; pass="cpass1_pass1"
  */
  ::Info("qaTrending::InitTrees::Begin","Detector %s, RefDet=%s",detector, referenceDet);
  AliExternalInfo info;
  Int_t treeCounter=0;
  // Load trees
  TObjArray * detArray=TString(referenceDet).Tokenize(";");
  Int_t nrefDets=detArray->GetEntries();
  TVectorF runCounter(5+nrefDets*2);  // <QADet>, <Logbook>, <Logbook.Det>, <MonAlisa>, <CPass>, <QA.RefDet>xnrefDets, <Logbook.RefDet>xnrefDets
  //
  ::Info("qaTrending::InitTrees::End","Laoding trees");
  treeQADet = info.GetTreeDataQA(detector,period, pass);
  if (!treeQADet){
    ::Error("qaTrending.InitTrees", "Input QA tree %s not available", detector); 
    return 0;
  }
  runCounter[treeCounter++]=treeQADet->Draw("run","1","goff");
  qaMap[TString::Format("QA.%s",detector).Data()]=treeQADet;
  //
  qaMap["Logbook"]=info.GetTree("Logbook",period,"");
  qaMap["Logbook"]->AddFriend(treeQADet,"QADet");
  treeQADet->AddFriend(qaMap["Logbook"],"Logbook");
  runCounter[treeCounter++]=treeQADet->Draw("run","1","goff");
  //
  qaMap["MonALISA.RCT"]=info.GetTree("MonALISA.RCT",period,pass);
  qaMap["MonALISA.RCT"]->AddFriend(treeQADet,"QADet");
  treeQADet->AddFriend(qaMap["MonALISA.RCT"],"MonALISA.RCT");
  runCounter[treeCounter++]=treeQADet->Draw("run","1","goff");
  //
  TTree *treeLogbookDetector =info.GetTree("Logbook.detector",period,"");
  if (treeLogbookDetector){
    if (detArray->GetEntries()>0){
      for (Int_t idet=0; idet<detArray->GetEntries(); idet++){
	// Load Logbook.RefDet
	const char *detName=detArray->At(idet)->GetName();	
	TTree * treeLog =treeLogbookDetector->CopyTree(TString::Format("detector==\"%s\"",detName).Data());
	if (treeLog->GetEntries()<=0){
	  ::Error("qaTrending.InitTrees","Missing Tree Logbook. %s - check the syntax",detName);
	}else{
	  treeLog->BuildIndex("run");
	  qaMap[TString::Format("Logbook.%s",detName).Data()]=treeLog;
	  treeLog->AddFriend(treeQADet, "QADet");
	  treeQADet->AddFriend(treeLog,  TString::Format("Logbook.%s",detName).Data());
	  runCounter[treeCounter++]=treeQADet->Draw("run","1","goff");
	}
	// Load QA.RefDet
	TTree * treeQARefDet = info.GetTreeDataQA(detName,period, pass);
	if (treeQARefDet){
	  qaMap[TString::Format("QA.%s",detName).Data()]=treeQARefDet;
	  treeQARefDet->AddFriend(treeQADet, "QADet");
	  treeQADet->AddFriend(treeQARefDet,  TString::Format("QA.%s",detName).Data());
	  runCounter[treeCounter++]=treeQADet->Draw("run","1","goff");
	}else{
	  ::Error("qaTrending.InitTrees","Missing Tree QA.%s - check the syntax",detName);
	}
      }
    }
  }
  //
  // Check consistency of data
  // 
  ::Info("qaTrending::InitTrees::End","Checking trees");
  TList *arrFriends = treeQADet->GetListOfFriends();
  for (Int_t ifriend=0; ifriend<arrFriends->GetEntries(); ifriend++){
    Int_t entries = treeQADet->Draw(TString::Format("run-%s.run", arrFriends->At(ifriend)->GetName()).Data(),"1","goff");
    Double_t mean=0;
    if (entries>0) {
      mean=TMath::Mean(entries, treeQADet->GetV1());
    }
    if (mean==0){
      ::Info("qaTrending::InitTrees", "Friend:\t%s\t%d\t%f", arrFriends->At(ifriend)->GetName(), entries,mean);
    }else{
      ::Error("qaTrending::InitTrees", "Friend:\t%s\t%d\t%f", arrFriends->At(ifriend)->GetName(), entries,mean);
    }
  }
  delete detArray;
  ::Info("qaTrending::InitTrees::End","Detector %s, RefDet=%s",detector, referenceDet);

  return treeQADet;  
}

void fit_and_weights_norm(){

    gROOT->ProcessLine(".L ~/cern/project/lhcbStyle.C");
    lhcbStyle();
    /*gStyle->SetLabelSize(0.05,"x");
    gStyle->SetLabelSize(0.05,"y");
    gStyle->SetTitleSize(0.05,"x");
    gStyle->SetPaperSize(20,26);
    gStyle->SetPadTopMargin(0.0);
    gStyle->SetPadRightMargin(0.05); // increase for colz plots
    gStyle->SetPadBottomMargin(0.0);
    gStyle->SetPadLeftMargin(0.14);
    gStyle->SetTitleH(0.01);*/
                                                                                    //
    const std::string filename("/afs/cern.ch/work/a/apmorris/private/cern/ntuples/new_tuples/normalisation_samples/Lb2JpsipK_2011_2012_signal_withbdt.root");           
    const std::string treename = "withbdt";                                         
    const std::string out_file_mass("~/cern/plots/fitting/Lb2JpsipK_2011_2012_mass_fit_after_bdtg3_cut_-055.pdf");                                   
                                                                                    //

    const std::string cuts("bdtg3>=-0.55");    
    
    TFile* file = TFile::Open( filename.c_str() );
    if( !file ) std::cout << "file " << filename << " does not exist" << std::endl;
    
    TTree* tree = (TTree*)file->Get( treename.c_str() );
    if( !tree ) std::cout << "tree " << treename << " does not exist" << std::endl;

    TTree* rTree1 = tree->CopyTree( cuts.c_str() );

    // -- signal, mass shape
    RooRealVar Lambda_b0_DTF_MASS_constr1("Lambda_b0_DTF_MASS_constr1","m(J/#psi p K^{-})", 5550., 5700., "MeV/c^{2}"); 
    RooRealVar Jpsi_M("Jpsi_M","m(#mu#mu)", 3000., 3200., "MeV/c^{2}"); 
    //RooRealVar chi_c_M("chi_c_M","m(J/#psi#gamma)", 3400., 3700., "MeV/c^{2}"); 
    RooRealVar mean("mean","mean", 5620., 5600., 5650.);
    RooRealVar sigma1("sigma1","sigma1", 10., 1., 100.);
    RooRealVar sigma2("sigma2","sigma2", 5.0, 1.0, 300.0);
    RooRealVar sigmaT("sigmaT", "sigmaT", 1.9, 1., 100.);
    RooRealVar alpha1("alpha1","alpha1", 1.0, 0.5, 5.0);
    RooRealVar n1("n1","n1", 1.5, 0.2, 15.0);
    RooRealVar alpha2("alpha2","alpha2", -0.5, -5.5, 0.0);
    RooRealVar n2("n2","n2", 1.5, 0.2, 10.0);
    RooRealVar nu("nu", "nu", 2., 0.7, 100.);
    //RooRealVar bkgcat_chic("bkgcat_chic","bkgcat_chic", 0, 100);
    RooRealVar bdtg3("bdtg3", "bdtg3", -1.0, 1.0);                                    //
    RooRealVar frac2("frac2","frac2", 0.3, 0., 1.);
    
    Lambda_b0_DTF_MASS_constr1.setBins(75);
    
    
    
    
    RooGaussian gauss1("gauss1","gauss1", Lambda_b0_DTF_MASS_constr1, mean, sigma1);
    RooGaussian gauss2("gauss2","gauss2", Lambda_b0_DTF_MASS_constr1, mean, sigma2);
    RooCBShape cb1("cb1","cb1", Lambda_b0_DTF_MASS_constr1, mean, sigma1, alpha1, n1); 
    RooCBShape cb2("cb2","cb2", Lambda_b0_DTF_MASS_constr1, mean, sigma2, alpha2, n2);
    
    //RooStudentT * student = new RooStudentT("student", "student", Lambda_b0_DTF_MASS_constr1, mean, sigmaT, nu);
     
    RooAddPdf sig("sig", "sig", RooArgList(cb1, cb2), RooArgList( frac2 ));
    RooRealVar cbRatio("cbRatio","cb Ratio", 0.3, 0.1, 1.0);
    RooRealVar sigYield("sigYield","sig Yield", 4e2, 1e1, 1e5);
    RooRealVar bgYield("bgYield","bg Yield", 1e2, 1e0, 5e5);

    //put in values from fit_MC here <<--- DON'T FORGET TO CHANGE THESE IF THE FIT CHANGES!!!
    /*
    EXT PARAMETER                                INTERNAL      INTERNAL
  NO.   NAME      VALUE            ERROR       STEP SIZE       VALUE
   1  alpha1       2.18020e+00   2.85078e-02   1.38432e-04  -2.56034e-01
   2  alpha2      -2.79102e+00   6.74385e-02   1.51818e-04  -1.49177e-02
   3  cbRatio      3.07172e-01   1.49204e-02   1.72642e-04  -5.69984e-01
   4  mean         5.61985e+03   9.58397e-03   5.56682e-05  -9.66293e-02
   5  n1           1.49358e+00   8.14447e-02   2.09300e-04  -9.70542e-01
   6  n2           1.45276e+00   1.09864e-01   2.59028e-04  -8.39538e-01
   7  sigma1       8.46303e+00   1.32851e-01   2.86985e-05  -1.01453e+00
   8  sigma2       4.93976e+00   3.42842e-02   5.03572e-06  -1.44512e+00
   
   
   
   
   4  mean         5.62097e+03   4.02152e-02   6.00497e-04   3.08772e-01
   5  sigYield     3.52933e+04   2.55400e+02   1.54032e-03  -1.69958e-02
   6  sigma1       1.22322e+01   1.10970e+00   2.87462e-03   1.63838e-01
   7  sigma2       5.54047e+00   1.41829e-01   1.08300e-03  -1.28653e-01
    */
    
   mean.setVal(5.62097e+03);
   sigma1.setVal(1.22322e+01);
   sigma2.setVal(5.54047e+00);
   
   mean.setConstant(true);
   sigma1.setConstant(true);
   sigma2.setConstant(true);
    
    
    
    //alpha1.setVal( 2.18020e+00 );
    //alpha2.setVal( -2.79102e+00 );
    //n1.setVal( 1.49358e+00 );
    //n2.setVal( 1.45276e+00 );
    //frac2.setVal( 3.81630e-01 );
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      probHistFi   = new TH1F( "MVA_Fisher_Proba",  "MVA_Fisher_Proba",  nbin, 0, 1 );
      rarityHistFi = new TH1F( "MVA_Fisher_Rarity", "MVA_Fisher_Rarity", nbin, 0, 1 );
   }

   // Prepare input tree (this must be replaced by your data source)
   // in this example, there is a toy tree with signal and one with background events
   // we'll later on use only the "signal" events for the test in this example.
   //   
   TFile *input(0);
   TString fname = "/tmp/chasco/ORIGINAL//Data_MuEG2011B_1.root";   
   if (!gSystem->AccessPathName( fname )) {
      input = TFile::Open( fname ); // check if file in local directory exists
   } 
   else { 
      input = TFile::Open( "http://root.cern.ch/files/tmva_class_example.root" ); // if not: download from ROOT server
   }
   
   if (!input) {
      std::cout << "ERROR: could not open data file" << std::endl;
      exit(1);
   }
   std::cout << "--- TMVAClassificationApp    : Using input file: " << input->GetName() << std::endl;
   
   //
   // prepare the tree
   // - here the variable names have to corresponds to your tree
   // - you can use the same variables as above which is slightly faster,
   //   but of course you can use different ones and copy the values inside the event loop
   //
   TTree* BigTree = (TTree*)input->Get("data");

TFile *tmp  = new TFile( "tmp.root","RECREATE" );

TTree* theTree = BigTree->CopyTree("((cat == 1) + (cat == 2))*(ln==0)*(Cosmic==0)*(fabs(Mass_Z - 91.18)<10)*(Pt_Z>30)*(DeltaPhi_metjet>0.5)*(Pt_J1 < 30)*(pfMEToverPt_Z > 0.4)*(pfMEToverPt_Z < 1.8)*((Pt_Jet_btag_CSV_max > 20)*(btag_CSV_max < 0.244) + (1-(Pt_Jet_btag_CSV_max > 20)))*(sqrt(pow(dilepPROJLong + 1.25*recoilPROJLong + 0.0*uncertPROJLong,2)*(dilepPROJLong + 1.25*recoilPROJLong + 0.0*uncertPROJLong > 0) + 1.0*pow(dilepPROJPerp + 1.25*recoilPROJPerp + 0.0*uncertPROJPerp,2)*(dilepPROJPerp + 1.25*recoilPROJPerp + 0.0*uncertPROJPerp > 0)) > 45.0)");
   std::cout << "--- Select signal sample" << std::endl;
   Float_t userVar1, userVar2;
//   theTree->SetBranchAddress( "var1", &userVar1 );
//   theTree->SetBranchAddress( "var2", &userVar2 );
//   theTree->SetBranchAddress( "var3", &var3 );
//   theTree->SetBranchAddress( "var4", &var4 );

theTree->SetBranchAddress( " Z_rapidity_z", &Z_rapidity_z);

theTree->SetBranchAddress( " THRUST_2D", &THRUST_2D);

theTree->SetBranchAddress( " L1_L2_cosangle", &L1_L2_cosangle);

theTree->SetBranchAddress( " TransMass_ZH150_uncl", &TransMass_ZH150_uncl);

theTree->SetBranchAddress( " TransMass_ZH150", &TransMass_ZH150);

theTree->SetBranchAddress( " DeltaPhi_ZH", &DeltaPhi_ZH);

theTree->SetBranchAddress( " DeltaPhi_ZH_uncl", &DeltaPhi_ZH_uncl);

theTree->SetBranchAddress( " CMAngle", &CMAngle);

theTree->SetBranchAddress( " CS_cosangle", &CS_cosangle);

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

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

   std::cout << "--- Processing: " << theTree->GetEntries() << " events" << std::endl;