UInt_t GetListOfTitles( TString & methodName, TList & titles, TDirectory *dir=0 )
   {
      // get the list of all titles for a given method
      // if the input dir is 0, gDirectory is used
      // returns a list of keys
      UInt_t ni=0;
      if (dir==0) dir = gDirectory;
      TDirectory* rfdir = (TDirectory*)dir->Get( methodName );
      if (rfdir==0) {
         cout << "+++ Could not locate directory '" << methodName << endl;
         return 0;
      }

      return GetListOfTitles( rfdir, titles );

      TList *keys = rfdir->GetListOfKeys();
      if (keys==0) {
         cout << "+++ Directory '" << methodName << "' contains no keys" << endl;
         return 0;
      }
      //
      TIter rfnext(rfdir->GetListOfKeys());
      TKey *rfkey;
      titles.Clear();
      titles.SetOwner(kFALSE);
      while ((rfkey = (TKey*)rfnext())) {
         // make sure, that we only look at histograms
         TClass *cl = gROOT->GetClass(rfkey->GetClassName());
         if (cl->InheritsFrom("TDirectory")) {
            titles.Add(rfkey);
            ni++;
         }
      }
      cout << "--- Found " << ni << " instance(s) of the method " << methodName << endl;
      return ni;
   }

void plotasym(int top){
  if (setup("vm","","pol__")==kFALSE) return;

  THStack* hs = new THStack(TString::Format("asym_top%d",top), TString::Format("asym_top%d",top));

  int itop = top-1;
  TFile* fy = _fyexp[itop];
  TIter nextkey(fy->GetListOfKeys());
  TKey *key;
  while (key = (TKey*)nextkey()) {
    TString Q2Wdirname = key->GetName();
    if(Q2Wdirname.EqualTo("hYW_Dir") || Q2Wdirname.EqualTo("hYW"))continue;
    cout << "Q2Wdirname = " << Q2Wdirname << endl;
    TString hname = TString::Format("%s/hAsym/Varset1/hAsym_ACC_CORR_phi",Q2Wdirname.Data());
    cout << "hname = " << hname << endl;
    TH1D* h = (TH1D*)fy->Get(hname);
    if (h==NULL) cout << "histogram not found" << endl;
    //h->Draw();
    hs->Add(h,"e1");
  }
  TCanvas *c = new TCanvas(hs->GetName(),hs->GetTitle());
  hs->Draw("pads");

}

void printBadEvents(string filename){
  string run = filename.substr(filename.find("_R000")+10, 1);
  int nrun = atoi(run.c_str());
  string fname(filename.substr(filename.find_last_of("/")+1) );
  cout << "FileName is " << fname.c_str() << endl;

  ofstream outfile;
  stringstream namefile;
  namefile << "BadEvents_" << nrun << ".txt";
  outfile.open(namefile.str().c_str());
  outfile << "Bad events in file " << fname.c_str() << endl; 

  TDirectory* topDir; 
  TFile* file = TFile::Open(filename.c_str());
  if (!file->IsOpen()) {
    cerr << "Failed to open " << filename << endl; 
    return;
  }
  string dir = "DQMData/Run " + run + "/ParticleFlow/Run summary/ElectronValidation/JetPtRes/BadEvents";
  topDir = dynamic_cast<TDirectory*>( file->Get(dir.c_str()));
  topDir->cd();
  if (topDir){
    TIter next(topDir->GetListOfKeys());
    TKey *key;
    while  ( (key = dynamic_cast<TKey*>(next())) ) {
      string sflag = key->GetName();
      string info(sflag.substr(1, sflag.find_first_of(">")-1 ) );
      string run(info.substr(0, info.find_first_of("_")) );
      string evt = info.substr( info.find_first_of("_")+1, info.find_last_of("_")-2);
      string ls(info.substr(info.find_last_of("_")+1,info.find_first_of("_")+1));
      string ptres = ( sflag.substr( sflag.find( "f=" ) + 2, 6 ) ).c_str();
      cout << "Event info: Run " << run << " LS " << ls << " Evt " << evt << " Jet Pt Res = " << ptres << endl; 
      outfile << "Event info: Run " << run << " LS " << ls << " Evt " << evt << " Jet Pt Res = " << ptres << endl; 
    }
  }
}

// input: - Input file (result from TMVA),
//        - use of TMVA plotting TStyle
void BDTControlPlots( TString fin = "TMVA.root", Bool_t useTMVAStyle = kTRUE )
{
   // set style and remove existing canvas'
   TMVAGlob::Initialize( useTMVAStyle );
  
   // checks if file with name "fin" is already open, and if not opens one
   TFile* file = TMVAGlob::OpenFile( fin );  

   // get all titles of the method BDT
   TList titles;
   UInt_t ninst = TMVAGlob::GetListOfTitles("Method_BDT",titles);
   if (ninst==0) {
      cout << "Could not locate directory 'Method_BDT' in file " << fin << endl;
      return;
   }
   // loop over all titles
   TIter keyIter(&titles);
   TDirectory *bdtdir;
   TKey *key;
   while ((key = TMVAGlob::NextKey(keyIter,"TDirectory"))) {
      bdtdir = (TDirectory *)key->ReadObj();
      bdtcontrolplots( bdtdir );
   }
}

 UInt_t GetListOfMethods( TList & methods, TDirectory *dir=0 )
 {
    // get a list of methods
    // the list contains TKey objects
    if (dir==0) dir = gDirectory;
    TIter mnext(dir->GetListOfKeys());
    TKey *mkey;
    methods.Clear();
    methods.SetOwner(kFALSE);
    UInt_t ni=0;
    while ((mkey = (TKey*)mnext())) {
       // make sure, that we only look at TDirectory with name Method_<xxx>
       TString name = mkey->GetClassName();
       TClass *cl = gROOT->GetClass(name);
       if (cl->InheritsFrom("TDirectory")) {
          if (TString(mkey->GetName()).BeginsWith("Method_")) {
             methods.Add(mkey);
             ni++;
          }
       }
    }
    cout << "--- Found " << ni << " classifier types" << endl;
    return ni;
 }

void drawBeamSpotGraphsAll (TFile* file, const char* fname=0)
{
  TDirectory* module = (TDirectory*)_file0->FindObjectAny("beamSpotFitPV");
  if ( module==0 )  return;

//   unsigned int nrun = module->GetListOfKeys()->GetSize();
//   TDirectory* currDir = gDirectory;
//   gROOT->cd();
//   TH1* hChi2 = new TH1F("runChi2","runChi2",nrun,0.5,nrun+0.5);
//   TH1* hConst = new TH1F("runConst","runConst",nrun,0.5,nrun+0.5);
//   TH1* hSlope = new TH1F("runSlope","runSlope",nrun,0.5,nrun+0.5);
//   currDir->cd();

  TIter iter(module->GetListOfKeys());
  TKey* key;
//   unsigned int irun(0);
//   float runSummary[6];
  while ( (key=(TKey*)iter()) ) {
    std::cout << key->GetName() << std::endl;
    drawBeamSpotGraphs(_file0,key->GetName(),fname);
    TIter ic(gROOT->GetListOfCanvases());
    TCanvas* c;
    while ( (c=(TCanvas*)ic()) )  delete c;
//     ++irun;
//     if ( runSummary[1]>0.1 )  hChi2->SetBinContent(irun,runSummary[0]/runSummary[1]);
//     hConst->SetBinContent(irun,runSummary[2]);
//     hConst->SetBinError(irun,runSummary[3]);
//     hSlope->SetBinContent(irun,runSummary[4]);
//     hSlope->SetBinError(irun,runSummary[5]);
//     hChi2->GetXaxis()->SetBinLabel(irun,key->GetName());
//     hConst->GetXaxis()->SetBinLabel(irun,key->GetName());
//     hSlope->GetXaxis()->SetBinLabel(irun,key->GetName());
//     if ( irun==3 )  return;
  }

}

    void ReadHist(TString currFName){
	    TFile *currFile = TFile::Open( currFName.Data(),"READONLY" );
    
        TIter next( gDirectory->GetListOfKeys() ); //-dir-content--

        Int_t oncePRsum=0;
        //--reading-flow-histograms-from-ROOT-file-------------(begin)----
        TKey *key; 
        while( (key = (TKey*)next()) ){ 
            //cout << " Classname " << key->GetClassName() << endl;
            //cout << " Name "      << key->GetName()      << endl; 
            //cout << " Title "     << key->GetTitle()     << endl;
            printf("Class=%s Name=%s Title=\"%s\"\n", key->GetClassName(), key->GetName(), key->GetTitle());

            TObject *obj = key->ReadObj();
            if ( obj->IsA()->InheritsFrom( TH1::Class() ) ) {
                TH1 *h = (TH1*)obj;
                //h = (TH1*)obj;
                //h->Draw();
                //c.Update();
                sscanf(key->GetTitle(), "%s %i %s %i %s %i"    , ch1, &iC, ch2, &iY, ch3, &iPt);
                //printf("%s - %i - %s - %i\n", ch1,  iC, ch2,  iY);
                TString s(key->GetName());
                if( s.Contains("hPhiPRp") ){ //hPhiPRp
                    printf("[PR]--%s[%i]%s[%i]%s[%i]--->\n",ch1,iC, ch2,iY, ch3,iPt);
                    hPhiPR[iC][iY][iPt] = new TH1F;
                    hPhiPR[iC][iY][iPt] = (TH1F*)obj;
                }
                if( s.Contains("hPhiABp") ){ //hPhiABp
                    printf("[AB]--%s[%i]%s[%i]%s[%i]--->\n",ch1,iC, ch2,iY, ch3,iPt);
                    hPhiAB[iC][iY][iPt] = new TH1F;
                    hPhiAB[iC][iY][iPt] = (TH1F*)obj;
                }
                if( s.Contains("hPhiEPp") ){ //hPhiEPp
                    printf("[EP]--%s[%i]%s[%i]%s[%i]--->\n",ch1,iC, ch2,iY, ch3,iPt);
                    hPhiEP[iC][iY][iPt] = new TH1F;
                    hPhiEP[iC][iY][iPt] = (TH1F*)obj;
                }
                if( s.Contains("hMulTrp") ){ //hMulTrp
                    printf("[Tr]--%s[%i]%s[%i]%s[%i]--->\n",ch1,iC, ch2,iY, ch3,iPt);
                    hMulTr[iC][iY][iPt] = new TH1F;
                    hMulTr[iC][iY][iPt] = (TH1F*)obj;
                }
            }
        }
    }

//________________________________________________________
void GFOverlay::Overlay(const TObjArray &dirs, const TObjArray &legends)
{
  // 'directories' must contain TDirectory and inheriting, being parallel with legends
  TDirectory *dir1 = 0;
  for (Int_t iDir = 0; !dir1 && iDir < dirs.GetEntriesFast(); ++iDir) {
    dir1 = static_cast<TDirectory*>(dirs[iDir]);
  }
  if (!dir1) return;

  const Int_t currentLayer = fLayer;
  fLayer += (fSummaries ? 2 : 1);
  std::vector<TH1*> meanHists, rmsHists;

  UInt_t counter = 0;
  TIter nextKey(dir1->GetListOfKeys());
//   //  while(TKey* key = static_cast <TKey*> (nextKey())) {
//   // OK, make CINT happy, i.e. make .L GFOverlay.C work without extending '+':
  TKey* key = NULL; 
  while ((key = static_cast <TKey*> (nextKey()))) {
    if (!fNames.IsEmpty() && !this->KeyContainsListMember(key->GetName(), fNames)) continue;
    if (this->KeyContainsListMember(key->GetName(), fSkipNames)) continue;

    TObjArray hists(this->GetTypeWithNameFromDirs(TH1::Class(), key->GetName(), dirs));
    if (this->AddHistsAt(hists, legends, currentLayer, counter) > 0) {
      if (fSummaries) {
	this->CreateFillMeanRms(hists, currentLayer, dir1->GetName(), meanHists, rmsHists);
      }
      ++counter;
    }
    
    TObjArray subDirs(this->GetTypeWithNameFromDirs(TDirectory::Class(), key->GetName(), dirs));
    if (subDirs.GetEntries()) { // NOT GetEntriesFast()!
      ::Info("GFOverlay::Overlay", "Key '%s' has directories to do recursion.", key->GetName());
      this->Overlay(subDirs, legends);
    }
  }
   
  // If mean/rms hists created, add them to manager:
  for (unsigned int iMean = 0; iMean < meanHists.size(); ++iMean) {
    fHistMan->AddHistSame(meanHists[iMean], currentLayer + 1, 0, legends[iMean]->GetName());
    fHistMan->AddHistSame(rmsHists[iMean], currentLayer + 1, 1, legends[iMean]->GetName());
  }
}

std::vector<std::string>
signalList(const char* dirName="", const char* pattern="", unsigned int debug=0)
{
  std::vector<std::string> histnames;
  TIter next(gDirectory->GetListOfKeys());
  TKey* iobj;
  unsigned int idx=0;
  while((iobj = (TKey*)next())){
    if(iobj->IsFolder()) continue;
    if(debug>2){ std::cout << "[" << ++idx << "] ...Found object: " << iobj->GetName() << " of type: " << iobj->GetClassName() << std::endl; }
    std::string fullpath(dirName); 
    fullpath += fullpath == std::string("") ? "" : "/"; fullpath += iobj->GetName();
    // why does \\w*_\\d+ not work to catch them all?!?
    if(std::string(pattern).empty() && httSignal(iobj->GetName())){
      histnames.push_back(fullpath);
    }
    else if(!std::string(pattern).empty() && match(iobj->GetName(), (char*)pattern)){
      histnames.push_back(fullpath);
    }
  }
  return histnames;
}

void MergeRootfile( char *core, TDirectory *target, TList *sourcelist ) {
  printf("merging only histos with core=%s=\n",core);

  //  cout << "Target path: " << target->GetPath() << endl;
  TString path( (char*)strstr( target->GetPath(), ":" ) );
  path.Remove( 0, 2 );

  TFile *first_source = (TFile*)sourcelist->First();
  first_source->cd( path );
  TDirectory *current_sourcedir = gDirectory;

  int nh=0;
  // loop over all keys in this directory
  TChain *globChain = 0;
  TIter nextkey( current_sourcedir->GetListOfKeys() );
  TKey *key;
  while ( (key = (TKey*)nextkey())) {
    const char *name=key->GetName();
    char * c1=strstr(name,core);
    if(c1==0) continue;
    if(c1-name>2) continue;
    nh++;
    if(nh%100==0) printf("nh=%d addingX %s\n",nh,name);
    
    // read object from first source file
    first_source->cd( path );
    TObject *obj = key->ReadObj();
    
    if ( obj->IsA()->InheritsFrom( "TH1" ) ) {
      // descendant of TH1 -> merge it
      
      //      cout << "Merging histogram " << obj->GetName() << endl;
      TH1 *h1 = (TH1*)obj;

      // loop over all source files and add the content of the
      // correspondant histogram to the one pointed to by "h1"
      TFile *nextsource = (TFile*)sourcelist->After( first_source );
      while ( nextsource ) {
        
        // make sure we are at the correct directory level by cd'ing to path
        nextsource->cd( path );
        TH1 *h2 = (TH1*)gDirectory->Get( h1->GetName() );
        if ( h2 ) {
          h1->Add( h2 );
          delete h2; // don't know if this is necessary, i.e. if 
                     // h2 is created by the call to gDirectory above.
        }

        nextsource = (TFile*)sourcelist->After( nextsource );
      }
    }
    else if ( obj->IsA()->InheritsFrom( "TTree" ) ) {
      
      // loop over all source files create a chain of Trees "globChain"
      const char* obj_name= obj->GetName();

      globChain = new TChain(obj_name);
      globChain->Add(first_source->GetName());
      TFile *nextsource = (TFile*)sourcelist->After( first_source );
      //      const char* file_name = nextsource->GetName();
      // cout << "file name  " << file_name << endl;
     while ( nextsource ) {
     	  
       globChain->Add(nextsource->GetName());
       nextsource = (TFile*)sourcelist->After( nextsource );
     }
     
    } else if ( obj->IsA()->InheritsFrom( "TDirectory" ) ) {
      // it's a subdirectory
      
      cout << "Found subdirectory " << obj->GetName() << endl;
      
      // create a new subdir of same name and title in the target file
      target->cd();
      TDirectory *newdir = target->mkdir( obj->GetName(), obj->GetTitle() );

      // newdir is now the starting point of another round of merging
      // newdir still knows its depth within the target file via
      // GetPath(), so we can still figure out where we are in the recursion
      MergeRootfile( core,newdir, sourcelist );
      
    } else {
      
      // object is of no type that we know or can handle
      cout << "Unknown object type, name: " 
           << obj->GetName() << " title: " << obj->GetTitle() << endl;
    }
    
    // now write the merged histogram (which is "in" obj) to the target file
    // note that this will just store obj in the current directory level,
    // which is not persistent until the complete directory itself is stored
    // by "target->Write()" below
    if ( obj ) {
      target->cd();
      
      //!!if the object is a tree, it is stored in globChain...
      if(obj->IsA()->InheritsFrom( "TTree" ))
	globChain->Write( key->GetName() );
      else
	obj->Write( key->GetName() );
//.........这里部分代码省略.........

void plot_efficiencies( TFile* file, Int_t type = 2, TDirectory* BinDir)
{
   // input:   - Input file (result from TMVA),
   //          - type = 1 --> plot efficiency(B) versus eff(S)
   //                 = 2 --> plot rejection (B) versus efficiency (S)

   Bool_t __PLOT_LOGO__  = kTRUE;
   Bool_t __SAVE_IMAGE__ = kTRUE;

   // the coordinates
   Float_t x1 = 0;
   Float_t x2 = 1;
   Float_t y1 = 0;
   Float_t y2 = 0.8;

   // reverse order if "rejection"
   if (type == 2) {
      Float_t z = y1;
      y1 = 1 - y2;
      y2 = 1 - z;    
      //      cout << "--- type==2: plot background rejection versus signal efficiency" << endl;
   }
   else {
      //  cout << "--- type==1: plot background efficiency versus signal efficiency" << endl;
   }
   // create canvas
   TCanvas* c = new TCanvas( "c", "the canvas", 200, 0, 650, 500 );

   // global style settings
   c->SetGrid();
   c->SetTicks();

   // legend
   Float_t x0L = 0.107,     y0H = 0.899;
   Float_t dxL = 0.457-x0L, dyH = 0.22;
   if (type == 2) {
      x0L = 0.15;
      y0H = 1 - y0H + dyH + 0.07;
   }
   TLegend *legend = new TLegend( x0L, y0H-dyH, x0L+dxL, y0H );
   legend->SetTextSize( 0.05 );
   legend->SetHeader( "MVA Method:" );
   legend->SetMargin( 0.4 );

   TString xtit = "Signal efficiency";
   TString ytit = "Background efficiency";  
   if (type == 2) ytit = "Background rejection";
   TString ftit = ytit + " versus " + xtit;

   if (TString(BinDir->GetName()).Contains("multicut")){
      ftit += "  Bin: ";
      ftit += (BinDir->GetTitle());
   }

   // draw empty frame
   if(gROOT->FindObject("frame")!=0) gROOT->FindObject("frame")->Delete();
   TH2F* frame = new TH2F( "frame", ftit, 500, x1, x2, 500, y1, y2 );
   frame->GetXaxis()->SetTitle( xtit );
   frame->GetYaxis()->SetTitle( ytit );
   TMVAGlob::SetFrameStyle( frame, 1.0 );

   frame->Draw();  

   Int_t color = 1;
   Int_t nmva  = 0;
   TKey *key, *hkey;

   TString hNameRef = "effBvsS";
   if (type == 2) hNameRef = "rejBvsS";

   TList hists;
   TList methods;
   UInt_t nm = TMVAGlob::GetListOfMethods( methods );
   //   TIter next(file->GetListOfKeys());
   TIter next(&methods);

   // loop over all methods
   while (key = (TKey*)next()) {
      TDirectory * mDir = (TDirectory*)key->ReadObj();
      TList titles;
      UInt_t ninst = TMVAGlob::GetListOfTitles(mDir,titles);
      TIter nextTitle(&titles);
      TKey *titkey;
      TDirectory *titDir;
      while ((titkey = TMVAGlob::NextKey(nextTitle,"TDirectory"))) {
         titDir = (TDirectory *)titkey->ReadObj();
         TString methodTitle;
         TMVAGlob::GetMethodTitle(methodTitle,titDir);
         TIter nextKey( titDir->GetListOfKeys() );
         while ((hkey = TMVAGlob::NextKey(nextKey,"TH1"))) {
            TH1 *h = (TH1*)hkey->ReadObj();    
            TString hname = h->GetName();
            if (hname.Contains( hNameRef ) && hname.BeginsWith( "MVA_" )) {
               h->SetLineWidth(3);
               h->SetLineColor(color);
               color++; if (color == 5 || color == 10 || color == 11) color++; 
               h->Draw("csame");
               hists.Add(h);
               nmva++;
            }
//.........这里部分代码省略.........

void compareplots_noRatio(){
  vector<TFile*> files; 
  files.push_back(new TFile("/storage/9/schweiger/analyseFxFx/pythia8/100kEvents/ttbarMergedFxFx8TeVCTEQ6M-extracted.root"));
  files.push_back(new TFile("/storage/9/schweiger/analyseFxFx/pythia8/100kEvents/mergingscale_30/ttbarMergedFxFxMS30GeVMECut10GeV8TeVCTEQ6M-extracted.root"));   
  files.push_back(new TFile("/storage/9/schweiger/analyseFxFx/pythia8/100kEvents/mergingscale_100/ttbarMergedMS100GeV8TeVCTEQ6M-extracted.root"));   
  files.push_back(new TFile("/storage/9/schweiger/analyseFxFx/pythia8/100kEvents/mergingscale_100/ttbarMergedMS100GeVMCCut50GeV8TeVCTEQ6M-extracted.root"));   
  


  
  vector<TString> names;
  names.push_back("ttbar +0/1 Jet, #mu_{Q}=10 GeV #mu_{ME} = 10 GeV");
  names.push_back("ttbar +0/1 Jet, #mu_{Q}=30 GeV #mu_{ME} = 10 GeV");
  names.push_back("ttbar +0/1 Jet, #mu_{Q}=100 GeV #mu_{ME} = 10 GeV");
  names.push_back("ttbar +0/1 Jet, #mu_{Q}=100 GeV #mu_{ME} = 50 GeV");
  
  
  vector<TString> titles;
  titles.push_back("Gen-Jet p_{T}  with pos weights (GeV)");
  titles.push_back("Gen-Jet p_{T} with neg weights (GeV)");
  titles.push_back("Gen-Jet p_{T} (GeV)");
  titles.push_back("Gen_Jet #phi with pos. weights");
  titles.push_back("Gen_Jet #phi with neg. weights");
  titles.push_back("Gen_Jet #phi");
  titles.push_back("Gen Jet #theta with pos weights");
  titles.push_back("Gen Jet #theta with neg weights");
  titles.push_back("Gen Jet #theta");
  titles.push_back("Gen Jet Energy with pos weights (GeV) ");
  titles.push_back("Gen Jet Energy with neg weights (GeV)");
  titles.push_back("Gen Jet Energy (GeV)");
  titles.push_back("p_{T} of hardest Gen-Jet with pos weights (GeV)");
  titles.push_back("p_{T} of hardest Gen-Jet with neg weights (GeV)");
  titles.push_back("p_{T} of hardest Gen-Jet (GeV)");
  titles.push_back("p_{T} of 2nd hardest Gen-Jet with pos weights (GeV)");
  titles.push_back("p_{T} of 2nd hardest Gen-Jet with neg weights (GeV)");
  titles.push_back("p_{T} of 2nd hardest Gen-Jet (GeV)");
  titles.push_back("#eta of hardest Gen-Jets with pos weights");
  titles.push_back("#eta of hardest Gen-Jets with neg weights");
  titles.push_back("#eta of hardest Gen-Jets");
  titles.push_back("Number of Gen-Jets with pos. weights");
  titles.push_back("Number of Gen-Jets with neg. weights");
  titles.push_back("Number of Gen-Jets");

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

// Show no statistics box
gStyle->SetOptStat(0);

TH1::SetDefaultSumw2();

// Main program part
  TIter nextkey(files.at(0)->GetListOfKeys());
  TKey *key;
  bool first=true;
  TCanvas* c = new TCanvas();
  c->Print("plots.pdf[");

  // Save also as pictures
  int pictureNumber = 0;

  int run = 0;
  while (key = (TKey*)nextkey()) {
    pictureNumber++;
    TString pictureName = TString::Format("%d.png",pictureNumber);


    vector<TH1F*> histos;
    histos.push_back((TH1F*)key->ReadObj());
    for(size_t i=1;i<files.size();i++){
      histos.push_back((TH1F*)files.at(i)->Get(histos.at(0)->GetName()));
    }
		       
    for(size_t i=0;i<histos.size();i++){
      if(i == 0){
	histos.at(i)->SetLineColor(kBlack);
      }
      if(i == 1){
	histos.at(i)->SetLineColor(kRed);
      }
      if(i == 2){
	histos.at(i)->SetLineColor(kBlue);
      }
      if(i == 3){
	histos.at(i)->SetLineColor(kGreen+2);
      }
      if(i == 4){
	histos.at(i)->SetLineColor(kMagenta-7);
      }
      if(i == 5){
	histos.at(i)->SetLineColor(kOrange+7);
      }
    }
   
    for(size_t i=0;i<histos.size();i++){
      histos.at(i)->Sumw2();
      histos.at(i)->Scale(1./histos.at(i)->Integral(),"width");
    }

// Set axis title
histos.at(0)->GetYaxis()->SetTitle("Normalized units"); 
//.........这里部分代码省略.........

void XMLimport::readKeyGroup( TKey * pParent )
{
    TKey * pT;
    if( pParent )
    {
        pT = new TKey( pParent, mpHost );
    }
    else
    {
        pT = new TKey( 0, mpHost );
    }
    mpHost->getKeyUnit()->registerKey( pT );
    pT->setIsActive( ( attributes().value("isActive") == "yes" ) );
    pT->mIsFolder = ( attributes().value("isFolder") == "yes" );
    if (module)
        pT->mModuleMember = true;

    while( ! atEnd() )
    {
        readNext();
        if( isEndElement() ) break;

        if( isStartElement() )
        {
            if( name() == "name" )
            {
                pT->mName = readElementText();
                continue;
            }
            else if( name() == "packageName")
            {
                pT->mPackageName = readElementText();
                continue;
            }
            else if( name() == "script")
            {
                QString script = readElementText();
                pT->setScript( script );
                continue;
            }
            else if( name() == "command")
            {
                pT->mCommand = readElementText();
                continue;
            }
            else if( name() == "keyCode" )
            {
                pT->mKeyCode = readElementText().toInt();
                continue;
            }
            else if( name() == "keyModifier" )
            {
                pT->mKeyModifier = readElementText().toInt();
                continue;
            }

            else if( name() == "KeyGroup" )
            {
                readKeyGroup( pT );
            }
            else if( name() == "Key" )
            {
                readKeyGroup( pT );
            }
            else
            {
                readUnknownKeyElement();
            }
        }
    }
}

// input: - Input file (result from TMVA)
//        - use of TMVA plotting TStyle
void mvas( TString fin = "TMVA.root", HistType htype = MVAType, Bool_t useTMVAStyle = kTRUE )
{
   // set style and remove existing canvas'
   TMVAGlob::Initialize( useTMVAStyle );

   // switches
   const Bool_t Save_Images     = kTRUE;

   // checks if file with name "fin" is already open, and if not opens one
   TFile* file = TMVAGlob::OpenFile( fin );  

   // define Canvas layout here!
   Int_t xPad = 1; // no of plots in x
   Int_t yPad = 1; // no of plots in y
   Int_t noPad = xPad * yPad ; 
   const Int_t width = 600;   // size of canvas

   // this defines how many canvases we need
   TCanvas *c = 0;

   // counter variables
   Int_t countCanvas = 0;

   // search for the right histograms in full list of keys
   TIter next(file->GetListOfKeys());
   TKey *key(0);   
   while ((key = (TKey*)next())) {

      if (!TString(key->GetName()).BeginsWith("Method_")) continue;
      if( ! gROOT->GetClass(key->GetClassName())->InheritsFrom("TDirectory") ) continue;

      TString methodName;
      TMVAGlob::GetMethodName(methodName,key);

      TDirectory* mDir = (TDirectory*)key->ReadObj();

      TIter keyIt(mDir->GetListOfKeys());
      TKey *titkey;
      while ((titkey = (TKey*)keyIt())) {
         if (!gROOT->GetClass(titkey->GetClassName())->InheritsFrom("TDirectory")) continue;

         TDirectory *titDir = (TDirectory *)titkey->ReadObj();
         TString methodTitle;
         TMVAGlob::GetMethodTitle(methodTitle,titDir);

         cout << "--- Found directory for method: " << methodName << "::" << methodTitle << flush;
         TString hname = "MVA_" + methodTitle;
         if      (htype == ProbaType  ) hname += "_Proba";
         else if (htype == RarityType ) hname += "_Rarity";
         TH1* sig = dynamic_cast<TH1*>(titDir->Get( hname + "_S" ));
         TH1* bgd = dynamic_cast<TH1*>(titDir->Get( hname + "_B" ));

         if (sig==0 || bgd==0) {
            if     (htype == MVAType)     
               cout << "mva distribution not available (this is normal for Cut classifier)" << endl;
            else if(htype == ProbaType)   
               cout << "probability distribution not available (this is normal for Cut classifier)" << endl;
            else if(htype == RarityType)  
               cout << "rarity distribution not available (this is normal for Cut classifier)" << endl;
            else if(htype == CompareType) 
               cout << "overtraining check not available (this is normal for Cut classifier)" << endl;
            else cout << endl;
         } 
         else {
            cout << endl;
            // chop off useless stuff
            sig->SetTitle( Form("TMVA response for classifier: %s", methodTitle.Data()) );
            if      (htype == ProbaType) 
               sig->SetTitle( Form("TMVA probability for classifier: %s", methodTitle.Data()) );
            else if (htype == RarityType) 
               sig->SetTitle( Form("TMVA Rarity for classifier: %s", methodTitle.Data()) );
            else if (htype == CompareType) 
               sig->SetTitle( Form("TMVA overtraining check for classifier: %s", methodTitle.Data()) );
         
            // create new canvas
            TString ctitle = ((htype == MVAType) ? 
                              Form("TMVA response %s",methodTitle.Data()) : 
                              (htype == ProbaType) ? 
                              Form("TMVA probability %s",methodTitle.Data()) :
                              (htype == CompareType) ? 
                              Form("TMVA comparison %s",methodTitle.Data()) :
                              Form("TMVA Rarity %s",methodTitle.Data()));
         
            TString cname = ((htype == MVAType) ? 
                             Form("output_%s",methodTitle.Data()) : 
                             (htype == ProbaType) ? 
                             Form("probability_%s",methodTitle.Data()) :
                             (htype == CompareType) ? 
                             Form("comparison_%s",methodTitle.Data()) :
                             Form("rarity_%s",methodTitle.Data()));

            c = new TCanvas( Form("canvas%d", countCanvas+1), ctitle, 
                             countCanvas*50+200, countCanvas*20, width, (Int_t)width*0.78 ); 
    
            // set the histogram style
            TMVAGlob::SetSignalAndBackgroundStyle( sig, bgd );
   
            // normalise both signal and background
//.........这里部分代码省略.........

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

    //  TFile *f = (*sourceList.begin())->getFile();
    TFile *f = (*itSourceFile)->getFile();
    //  TTree *tree= (*sourceList.begin())->getTree();
    TTree *tree= (*itSourceFile)->getTree();
    int myLineColor = (*itSourceFile)->getLineColor();
    int myLineStyle = (*itSourceFile)->getLineStyle();
    TString myLegendName = (*itSourceFile)->getName();
    if (!f || !tree) {
      std::cout << "PlotAlignmentValidation::addHists: no tree or no file" << std::endl;
      return 0;
    }

    // Todo: TLegend?
  
    // first loop on tree to find out which entries (i.e. modules) fulfill the selection
    // 'Entry$' gives the entry number in the tree
    Long64_t nSel = tree->Draw("Entry$", selection, "goff");
    if (nSel == -1) return 0; // error in selection
    if (nSel == 0) {
      std::cout << "PlotAlignmentValidation::addHists: no selected module." << std::endl;
      return 0;
    }
    // copy entry numbers that fulfil the selection
    const std::vector<double> selected(tree->GetV1(), tree->GetV1() + nSel);

    TH1 *h = 0;       // becomes result
    UInt_t nEmpty = 0;// selected, but empty hists
    Long64_t nentries =  tree->GetEntriesFast();
    std::vector<double>::const_iterator iterEnt = selected.begin();

    // second loop on tree:
    // for each selected entry get the hist from the file and merge
    TkOffTreeVariables *treeMem = 0; // ROOT will initialise
    tree->SetBranchAddress("TkOffTreeVariables", &treeMem);
    for (Long64_t i = 0; i < nentries; i++){
      if (i < *iterEnt - 0.1             // smaller index (with tolerance): skip
	  || iterEnt == selected.end()) { // at the end: skip 
	continue;
      } else if (TMath::Abs(i - *iterEnt) < 0.11) {
	++iterEnt; // take this entry!
      } else std::cout << "Must not happen: " << i << " " << *iterEnt << std::endl;

      tree->GetEntry(i);
      if (printModuleIds) {
	std::cout << treeMem->moduleId << ": " << treeMem->entries << " entries" << std::endl;
      }
      if (treeMem->entries <= 0) {  // little speed up: skip empty hists
	++nEmpty;
	continue;
      }
      TString hName;
      switch(rType) {
      case xPrimeRes:     hName = treeMem->histNameX.c_str();          break;
      case yPrimeRes:     hName = treeMem->histNameY.c_str();          break;
      case xPrimeNormRes: hName = treeMem->histNameNormX.c_str();      break;
      case yPrimeNormRes: hName = treeMem->histNameNormY.c_str();      break;
      case xRes:          hName = treeMem->histNameLocalX.c_str();     break;
      case yRes:          hName = treeMem->histNameLocalY.c_str();     break;
      case xNormRes:      hName = treeMem->histNameNormLocalX.c_str(); break;
	/*case yResNorm:      hName = treeMem->histNameNormLocalY.c_str(); break;*/
      case ResXvsXProfile: hName = treeMem->profileNameResXvsX.c_str();    break;
      case ResXvsYProfile: hName = treeMem->profileNameResXvsY.c_str();    break;
      case ResYvsXProfile: hName = treeMem->profileNameResYvsX.c_str();    break;
      case ResYvsYProfile: hName = treeMem->profileNameResYvsY.c_str();    break;
      }
      TKey *histKey = f->FindKeyAny(hName);
      TH1 *newHist = (histKey ? static_cast<TH1*>(histKey->ReadObj()) : 0);
      if (!newHist) {
	std::cout << "Hist " << hName << " not found in file, break loop." << std::endl;
	break;
      }
      newHist->SetLineColor(myLineColor);
      newHist->SetLineStyle(myLineStyle);
      if (!h) { // first hist: clone, but rename keeping only first part of name
	TString name(newHist->GetName());
	Ssiz_t pos_ = 0;
	for (UInt_t i2 = 0; i2 < 3; ++i2) pos_ = name.Index("_", pos_+1);
	name = name(0, pos_); // only up to three '_'
	h = static_cast<TH1*>(newHist->Clone("summed_"+name));
	//      TString myTitle = Form("%s: %lld modules", selection, nSel);
	//	h->SetTitle( myTitle );
      } else { // otherwise just add
	h->Add(newHist);
      }
      delete newHist;
    }

    std::cout << "PlotAlignmentValidation::addHists" << "Result is merged from " << nSel-nEmpty
	      << " modules, " << nEmpty << " hists were empty." << std::endl;

    if (nSel-nEmpty == 0) continue;

    myLegend->AddEntry(myLegendName, myLegendName, "L");
    
    retHistoStack->Add(h);
  }
  myLegend->Draw();
  return retHistoStack;
}

void th2polyEurope(Int_t npoints=500000)
{
   Int_t i,j;
   Double_t lon1 = -25;
   Double_t lon2 =  35;
   Double_t lat1 =  34;
   Double_t lat2 =  72;
   Double_t R = (lat2-lat1)/(lon2-lon1);
   Int_t W    = 800;
   Int_t H    = (Int_t)(R*800);
   gStyle->SetStatX(0.28);
   gStyle->SetStatY(0.45);
   gStyle->SetStatW(0.15);

   // Canvas used to draw TH2Poly (the map)
   TCanvas *ce = new TCanvas("ce", "ce",0,0,W,H);
   ce->ToggleEventStatus();
   ce->SetGridx();
   ce->SetGridy();

   // Real surfaces taken from Wikipedia.
   const Int_t nx = 36;
   // see http://en.wikipedia.org/wiki/Area_and_population_of_European_countries
   const char *countries[nx] = {
      "france",     "spain",  "sweden",  "germany",       "finland",
      "norway",     "poland", "italy",   "yugoslavia",    "united_kingdom",
      "romania",    "belarus","greece",  "czechoslovakia","bulgaria",
      "iceland",    "hungary","portugal","austria",       "ireland",
      "lithuania",  "latvia", "estonia", "denmark",       "netherlands",
      "switzerland","moldova","belgium", "albania",       "cyprus",
      "luxembourg", "andorra","malta",   "liechtenstein", "san_marino", "monaco" };
   Float_t surfaces[nx] = {
      547030,        505580,   449964,      357021,        338145,
      324220,        312685,   301230,      255438,        244820,
      237500,        207600,   131940,      127711,        110910,
      103000,         93030,    89242,       83870,         70280,
      65200,         64589,    45226,       43094,         41526,
      41290,         33843,    30528,       28748,          9250,
      2586,           468,      316,         160,            61, 2};

   TH1F *h = new TH1F("h","Countries surfaces (in km^{2})",3,0,3);
   for (i=0; i<nx; i++) h->Fill(countries[i], surfaces[i]);
   h->LabelsDeflate();

   TFile::SetCacheFileDir(".");
   TFile *f;
   f = TFile::Open("http://root.cern.ch/files/europe.root","cacheread");

   if (!f) {
      printf("Cannot access europe.root. Is internet working ?\n");
      return;
   }

   TH2Poly *p = new TH2Poly(
             "Europe",
             "Europe (bin contents are normalized to the surfaces in km^{2})",
             lon1,lon2,lat1,lat2);
   p->GetXaxis()->SetNdivisions(520);
   p->GetXaxis()->SetTitle("longitude");
   p->GetYaxis()->SetTitle("latitude");

   p->SetContour(100);

   TMultiGraph *mg;
   TKey *key;
   TIter nextkey(gDirectory->GetListOfKeys());
   while ((key = (TKey*)nextkey())) {
      TObject *obj = key->ReadObj();
      if (obj->InheritsFrom("TMultiGraph")) {
         mg = (TMultiGraph*)obj;
         p->AddBin(mg);
      }
   }

   TRandom r;
   Double_t longitude, latitude;
   Double_t x, y, pi4 = TMath::Pi()/4, alpha = TMath::Pi()/360;

   gBenchmark->Start("Partitioning");
   p->ChangePartition(100, 100);
   gBenchmark->Show("Partitioning");

   // Fill TH2Poly according to a Mercator projection.
   gBenchmark->Start("Filling");
   for (i=0; i<npoints; i++) {
      longitude = r.Uniform(lon1,lon2);
      latitude  = r.Uniform(lat1,lat2);
      x         = longitude;
      y         = 38*TMath::Log(TMath::Tan(pi4+alpha*latitude));
      p->Fill(x,y);
   }
   gBenchmark->Show("Filling");

   Int_t nbins = p->GetNumberOfBins();
   Double_t maximum = p->GetMaximum();


   // h2 contains the surfaces computed from TH2Poly.
   TH1F *h2 = (TH1F *)h->Clone("h2");
   h2->Reset();
//.........这里部分代码省略.........

void plotauto(TString infil