TTree *ToyTree(TString dirname="test/dato/fitres/Hgg_Et-toys/0.01-0.00", TString fname="outProfile-scaleStep2smearing_7-Et_25-trigger-noPF-EB.root", TString opt="", int nSmooth=10){
  TString outDir=dirname; outDir.ReplaceAll("fitres","img");
  outDir="tmp/";
  //std::map<TString, TH2F *> deltaNLL_map;

  //bool smooth=false;
  //if(opt.Contains("smooth")) smooth=true;
  

  /*------------------------------ Plotto */
  TCanvas c("ctoy","c");
  
  
  TTree *toys = new TTree("toys","");
  toys->SetDirectory(0);
  Double_t constTerm_tree, constTermTrue_tree;
  Double_t alpha_tree, alphaTrue_tree;
  char catName[100]; 
  Int_t catIndex;
  toys->Branch("constTerm", &constTerm_tree, "constTerm/D");
  toys->Branch("alpha", &alpha_tree, "alpha/D");
  toys->Branch("constTermTrue", &constTermTrue_tree, "constTermTrue/D");
  toys->Branch("alphaTrue", &alphaTrue_tree, "alphaTrue/D");

  toys->Branch("catName", catName, "catName/C");
  toys->Branch("catIndex", &catIndex, "catIndex/I");
  std::map<TString, Int_t> catIndexMap;

  ///1/
  for(int itoy =2; itoy <= 50; itoy++){
    TString filename=dirname+"/"; filename+=itoy; filename+="/"+fname;
    TString fout=dirname+"/"; fout+=itoy; fout+="/";    
    TFile f_in(filename, "read");
    if(f_in.IsZombie()){
      std::cerr << "File opening error: " << filename << std::endl;
      continue; //return NULL;
    }
    //std::cout << filename << std::endl;
    TList *KeyList = f_in.GetListOfKeys();
    //std::cout << KeyList->GetEntries() << std::endl;
    for(int i =0; i <  KeyList->GetEntries(); i++){
      c.Clear();
      TKey *key = (TKey *)KeyList->At(i);
      if(TString(key->GetClassName())!="RooDataSet") continue;
      RooDataSet *dataset = (RooDataSet *) key->ReadObj();
    
      TString constTermName = dataset->GetName();
      TString alphaName=constTermName; alphaName.ReplaceAll("constTerm","alpha");
      if(constTermName.Contains("scale")) continue;
      if(constTermName.Contains("alpha")) continue;
      if(constTermName.Contains("1.4442-gold")) continue;
      TTree *tree = dataset2tree(dataset);

      TGraph *rhoGraph = GetRho(tree, alphaName, constTermName);

      rhoGraph->SaveAs(fout+"rhoGraph-"+constTermName+".root");


      TGraphErrors bestFit_ = bestFit(tree, alphaName, constTermName);
      //TString binning="(241,-0.0005,0.2405,61,-0.00025,0.03025)"; //"(40,0.00025,0.02025,61,-0.0022975,0.1401475)";
      TString binning="(241,-0.0005,0.2405,301,-0.00005,0.03005)"; 

      TH2F *hist = prof2d(tree, constTermName, alphaName, "nll", binning, true, nSmooth, opt);
      //hist->SaveAs("myhist.root");
      
      Int_t iBinX, iBinY;
      hist->GetBinWithContent2(0.0002,iBinX,iBinY,1,-1,1,-1,0.0000001);
	
      //      if(iBinX!=0 && iBinY!=0 && iBinX < 41 && iBinY < 62){
      {
	TString catName_=constTermName; catName_.ReplaceAll("constTerm_",""); catName_.ReplaceAll("-","_");
	if(catIndexMap.count(catName_)==0) catIndexMap.insert(std::pair<TString,Int_t>(catName_,catIndexMap.size()));
	catIndex=catIndexMap[catName_];	
	constTerm_tree =  hist->GetYaxis()->GetBinCenter(iBinY);
	alpha_tree = hist->GetXaxis()->GetBinCenter(iBinX);
	sprintf(catName,"%s", catName_.Data());
	bestFit_.GetPoint(0, constTermTrue_tree,alphaTrue_tree);
// 	std::cout << constTerm_tree << " " << constTermTrue_tree 
// 		  << "\t" << alpha_tree << " " << alphaTrue_tree 
// 		  << std::endl;

	if(opt.Contains("scandiff")){
	  constTermTrue_tree = getMinimumFromTree(tree, "nll",TString(constTermName).ReplaceAll("-","_"));
	} else       if(opt.Contains("scan")){
	  constTerm_tree = getMinimumFromTree(tree, "nll",TString(constTermName).ReplaceAll("-","_"));
	}
	//std::cout << iBinX << "\t" << iBinY << "\t" << constTerm_tree - getMinimumFromTree(tree, "nll",TString(constTermName).ReplaceAll("-","_")) << std::endl;
	
	toys->Fill();
//       }else{
// 	hist->SaveAs("myhist.root");
// 	exit(0);
      }

      
      delete tree;
      delete hist;
    
    }
    f_in.Close();
//.........这里部分代码省略.........

TProfile *extractFlowVZEROsingle(Int_t icentr,Int_t spec,Int_t arm,Bool_t isMC,Float_t pTh,Int_t addbin,const char *nameSp,Float_t detMin,Float_t detMax,Int_t chMin,Int_t chMax){
  LoadLib();

  pTh += 0.00001;

  // NUA correction currently are missing
  char name[100];
  char stringa[200];

  snprintf(name,100,"AnalysisResults.root");
  if(!fo) fo = new TFile(name);
  snprintf(name,100,"contVZEROv%i",arm);
  TList *cont = (TList *) fo->Get(name);

  cont->ls();

  Float_t xMin[5] = {icentr/*centrality bin*/,chMin/*charge*/,pTh/*prob*/,-TMath::Pi()/arm/*Psi*/,detMin/*PID mask*/};
  Float_t xMax[5] = {icentr+addbin,chMax,1.0,TMath::Pi()/arm,detMax};

  cont->ls();

  TProfile *p1 = cont->At(2);
  TProfile *p2 = cont->At(3);
  TProfile *p3 = cont->At(4);

  TH2F *hPsi2DA = cont->At(5);
  TH2F *hPsi2DC = cont->At(6);

  TH1D *hPsiA = hPsi2DA->ProjectionY("PsiA",icentr+1,icentr+addbin+1);
  TH1D *hPsiC = hPsi2DC->ProjectionY("PsiC",icentr+1,icentr+addbin+1);

  if(!fPsi) fPsi = new TF1("fPsi","pol0",-TMath::Pi()/arm,TMath::Pi()/arm);
  hPsiA->Fit(fPsi,"0");
  Float_t offsetA = fPsi->GetParameter(0);
  hPsiC->Fit(fPsi,"0");
  Float_t offsetC = fPsi->GetParameter(0);

  Int_t nbinPsi = hPsiA->GetNbinsX();

  Float_t *NUAcorrA = new Float_t[nbinPsi];
  Float_t *NUAcorrC = new Float_t[nbinPsi];

  for(Int_t i=0;i < nbinPsi;i++){
    NUAcorrA[i] = offsetA/(hPsiA->GetBinContent(i+1));
    NUAcorrC[i] = offsetC/(hPsiC->GetBinContent(i+1));
  }

  Float_t res1=0,res2=0,res3=0; 
  Float_t eres1=0,eres2=0,eres3=0; 

  for(Int_t i = icentr; i <= icentr+addbin;i++){
    if(p1->GetBinError(i+1)){
      eres1 += 1./p1->GetBinError(i+1)/p1->GetBinError(i+1);
      res1 += p1->GetBinContent(i+1)/p1->GetBinError(i+1)/p1->GetBinError(i+1);      
    }
    if(p2->GetBinError(i+1)){
      eres2 += 1./p2->GetBinError(i+1)/p2->GetBinError(i+1);
      res2 += p2->GetBinContent(i+1)/p2->GetBinError(i+1)/p2->GetBinError(i+1);      
    }
    if(p3->GetBinError(i+1)){
      eres3 += 1./p3->GetBinError(i+1)/p3->GetBinError(i+1);
      res3 += p3->GetBinContent(i+1)/p3->GetBinError(i+1)/p3->GetBinError(i+1);      
    }
  }

  res1 /= eres1;
  res2 /= eres2;
  res3 /= eres3;
  
  eres1 = sqrt(1./eres1);
  eres2 = sqrt(1./eres2);
  eres3 = sqrt(1./eres3);

  AliFlowVZEROResults *a = (AliFlowVZEROResults *) cont->At(0);
  AliFlowVZEROResults *b = (AliFlowVZEROResults *) cont->At(1);
  TProfile *pp,*pp2;
  if(kNUAcorr){ // with NUA corrections
    pp = a->GetV2reweight(spec,xMin,xMax,3,NUAcorrA);
    pp2 = b->GetV2reweight(spec,xMin,xMax,3,NUAcorrC);
  }
  else{
    pp = a->GetV2(spec,xMin,xMax);
    pp2 = b->GetV2(spec,xMin,xMax);
  }
  
  Float_t scaling = sqrt(res1*res3/res2);
  if(kVZEROrescorr){
    pp->Scale(1./scaling);
  }

  Float_t err1_2 = eres1*eres1/res1/res1/4 +
    eres2*eres2/res2/res2/4 +
    eres3*eres3/res3/res3/4;
  Float_t err2_2 = err1_2;
  err1_2 /= scaling*scaling;
  printf("resolution V0A = %f +/- %f\n",scaling,err1_2);
  scaling = sqrt(res2*res3/res1);
  err2_2 /= scaling*scaling;
  if(kVZEROrescorr){
    pp2->Scale(1./scaling);
//.........这里部分代码省略.........

void ContourList(){
 
   const Double_t PI = TMath::Pi(); 
    
   TCanvas* c = new TCanvas("c","Contour List",0,0,600,600);
   c->SetRightMargin(0.15);
   c->SetTopMargin(0.15);
 
   Int_t i, j, TotalConts;
    
   Int_t nZsamples   = 80;
   Int_t nPhiSamples = 80;

   Double_t HofZwavelength = 4.0;       // 4 meters 
   Double_t dZ             =  HofZwavelength/(Double_t)(nZsamples - 1);
   Double_t dPhi           = 2*PI/(Double_t)(nPhiSamples - 1);
 
   TArrayD z(nZsamples);
   TArrayD HofZ(nZsamples);
   TArrayD phi(nPhiSamples);
   TArrayD FofPhi(nPhiSamples);

   
   // Discretized Z and Phi Values
   for ( i = 0; i < nZsamples; i++) {
      z[i] = (i)*dZ - HofZwavelength/2.0;
      HofZ[i] = SawTooth(z[i], HofZwavelength);
   }

   for(Int_t i=0; i < nPhiSamples; i++){
      phi[i] = (i)*dPhi;      
      FofPhi[i] = sin(phi[i]); 
   }

   // Create Histogram
   TH2D *HistStreamFn = new TH2D("HstreamFn", 
   "#splitline{Histogram with negative and positive contents. Six contours are defined.}{It is plotted with options CONT LIST to retrieve the contours points in TGraphs}", 
   nZsamples, z[0], z[nZsamples-1], nPhiSamples, phi[0], phi[nPhiSamples-1]);
     
   // Load Histogram Data
   for (Int_t i = 0; i < nZsamples; i++) {
      for(Int_t j = 0; j < nPhiSamples; j++){
         HistStreamFn->SetBinContent(i,j, HofZ[i]*FofPhi[j]);
      } 
   }

   gStyle->SetPalette(1); 
   gStyle->SetOptStat(0);
   gStyle->SetTitleW(0.99);
   gStyle->SetTitleH(0.08);
   
   Double_t contours[6];
   contours[0] = -0.7;
   contours[1] = -0.5;
   contours[2] = -0.1;
   contours[3] =  0.1;
   contours[4] =  0.4;
   contours[5] =  0.8;
     
   HistStreamFn->SetContour(6, contours); 
 
   // Draw contours as filled regions, and Save points
   HistStreamFn->Draw("CONT Z LIST");
   c->Update(); // Needed to force the plotting and retrieve the contours in TGraphs
    
   // Get Contours
   TObjArray *conts = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
   TList* contLevel = NULL;
   TGraph* curv     = NULL;    
   Int_t nGraphs    = 0;
   Int_t TotalConts = 0;
    
   
   if (conts == NULL){
      printf("*** No Contours Were Extracted!\n");
      TotalConts = 0;
      return;
   } else {
      TotalConts = conts->GetSize();
   }

   printf("TotalConts = %d\n", TotalConts);
        
   for(i = 0; i < TotalConts; i++){       
      contLevel = (TList*)conts->At(i);
      printf("Contour %d has %d Graphs\n", i, contLevel->GetSize());
      nGraphs += contLevel->GetSize();
   }
    
   nGraphs = 0;
    
   TCanvas* c1 = new TCanvas("c1","Contour List",610,0,600,600);
   c1->SetTopMargin(0.15);
   TH2F *hr = new TH2F("hr",
   "#splitline{Negative contours are returned first (highest to lowest). Positive contours are returned from}{lowest to highest. On this plot Negative contours are drawn in red and positive contours in blue.}",
   2, -2, 2, 2, 0, 6.5);

   hr->Draw();
   Double_t x0, y0, z0;
   TLatex l;
//.........这里部分代码省略.........

void makePlots( const char * inputA, const char * inputB, const char * option) 
{
  
  //Output path
  TString path("./paper01-plots/probs/");
  
  TString dataPxxA = TString( "EarthA" ) + TString("_") + TString( "0" ) + TString("_") + TString("Pme") +  TString("/data");
  TString dataPxxB = TString( "EarthB" ) + TString("_") + TString( "0" ) + TString("_") + TString("Pme") +  TString("/data");
  
  TList * v_Labels = new TList();
  TObjString *label;
  label = new TObjString( "Earth: Step function" );
  v_Labels->Add( label ); 
  label = new TObjString( "Earth: Realistic" );
  v_Labels->Add( label ); 
  
  TFile * f1 = new TFile(inputA);
  TFile * f2 = new TFile(inputB);
  
  f1->cd();
  TTree * PxxATreeNu = (TTree*)gDirectory->Get( dataPxxA.Data() );
    
  f2->cd();
  TTree * PxxBTreeNu = (TTree*)gDirectory->Get( dataPxxB.Data() );
  
  //Branches
  double xx = 0.0;
  double yy = 0.0;
  
  TString cname = TString("Earth") + TString("_") + TString("Pee");
  
  TGraph * ProbNu[3];
  ProbNu[0] = new TGraph();
  ProbNu[1] = new TGraph();
  
  PxxATreeNu->SetBranchAddress("Ex",&xx);
  PxxATreeNu->SetBranchAddress("Pb",&yy);
  
  Long64_t nentries = PxxATreeNu->GetEntries();
  
  for (Long64_t i=0;i<nentries;i++) {
    PxxATreeNu->GetEntry(i);
    ProbNu[0]->SetPoint( i, xx, yy);
    if( i < 5 ) 
      std::cout << i << " " << xx << " " << yy << std::endl;
  }
  
  PxxBTreeNu->SetBranchAddress("Ex",&xx);
  PxxBTreeNu->SetBranchAddress("Pb",&yy);
  
  nentries = PxxBTreeNu->GetEntries();
  
  for (Long64_t i=0;i<nentries;i++) {
    PxxBTreeNu->GetEntry(i);
    ProbNu[1]->SetPoint( i, xx, yy);
  }
  
  int ndataset = 0;
  TList * allgraphs = new TList(); //this is a ROOT container. It will store all of your Graphs
  
  allgraphs->Add( ProbNu[0] ); 
  ndataset++;

  allgraphs->Add( ProbNu[1] ); 
  ndataset++;
  
  //Datasets options (Markers: style, color, size) : You can also do it by hand using the interactive Editor
  int style[5];
  int color[5];
  float size[5];
  
  //For dataset No1
  style[0] = 21;
  color[0] = 2;
  size[0]  = 0.7;
  
  //For dataset No2
  style[1] = 20;
  color[1] = 4;
  size[1]  = 0.8;
  
  //For dataset No3
  style[2] = 22;
  color[2] = 38;
  size[2]  = 0.8;

  TCanvas * c1 = new TCanvas( cname.Data(), "Oscillation probabilities", 184,60,861,263); 
  c1->cd();
  gPad->SetGridx();
  gPad->SetGridy();
  gPad->SetLogx();

  float ymin = 0.0;
  float ymax = 0.5;
  float xmin = 0.9e9;
  float xmax = 1.0e10;

  TLegend * leg = new TLegend(0.14,0.58,0.35,0.85);
    
  leg->SetBorderSize(0);
//.........这里部分代码省略.........

void MakePIDqaReport(const char* inputFile, const char* outputFile="PIDqaReport.pdf", TString dirInFile = "")
{
  //
  // Make a pdf file with the efficiency report
  //

  LoadLibs();
  SetupStyle();

  TFile f(inputFile);
  if (!f.IsOpen()){
    printf("Could not open file '%s'\n",f.GetName());
    return;
  }
  
  TString listName = "PIDqa";
  if (dirInFile != "")
    listName = listName.Prepend(Form("%s/", dirInFile.Data()));
  
  printf("%s", listName.Data());
  TList *qaList = (TList*) f.Get(listName.Data());
  if (!qaList){
    printf("Could not find list '%s' in file '%s'\n",listName.Data(), f.GetName());
    return;
  }

  fCanvas=new TCanvas;

  TPDF p(outputFile);

  //
  // Invariant mass plots
  //


  //
  // Make QA info
  //

  // ITS PID
  PublishCanvas(qaList,"ITS","hNsigmaP_ITS_%s");

  // TPC PID
  TList *qaListTPC = (TList*)qaList->FindObject("TPC");
  if (qaListTPC){
    PublishCanvas(qaListTPC,"TPCBasic","hNsigmaP_TPC_Basic_%s");
    PublishCanvas(qaListTPC,"TPCV0","hNsigmaP_TPC_V0_%s");
    //   if (man->GetCurrentPeriod()=="11h"){
      //     PublishCanvas(qaListTPC,"TPC","hNsigmaP_TPC_Basic_%s_Hybrid","Hybrid");
    //     PublishCanvas(qaListTPC,"TPC","hNsigmaP_TPC_Basic_%s_OROChigh","OROChigh");
    //   }
  }
  else {
    printf("Could not find list '%s/TPC' in file '%s'\n", listName.Data(), f.GetName());
  }

  // TPC PID after 3 sigma TOF cut
  PublishCanvas(qaList,"TPC_TOF","hNsigmaP_TPC_TOF_%s");

  // TOF PID
  PublishCanvas(qaList,"TOF","hNsigmaP_TOF_%s");

  // TRD PID
  fCanvas->Divide(2,3);
  TH2 *hLikeP_TRD_3tls_electron=Get2DHistogramfromList(qaList,"TRD","hLikeP_TRD_3tls_electron");
  TH2 *hLikeP_TRD_3tls_pion=Get2DHistogramfromList(qaList,"TRD","hLikeP_TRD_3tls_pion");
  TH2 *hLikeP_TRD_4tls_electron=Get2DHistogramfromList(qaList,"TRD","hLikeP_TRD_4tls_electron");
  TH2 *hLikeP_TRD_4tls_pion=Get2DHistogramfromList(qaList,"TRD","hLikeP_TRD_4tls_pion");
  TH2 *hLikeP_TRD_5tls_electron=Get2DHistogramfromList(qaList,"TRD","hLikeP_TRD_5tls_electron");
  TH2 *hLikeP_TRD_5tls_pion=Get2DHistogramfromList(qaList,"TRD","hLikeP_TRD_5tls_pion");

  /*
   *  cTRDnsigma[countcanvas]->cd(1);
   *  TPaveText pt3TRD(.02,.02,.49,.52);
   *  pt3TRD.SetTextAlign(11);
   *  pt3TRD.SetTextSizePixels(16);
   *  pt3TRD.AddText(Form(" TRD PID QA %s.%s.%d", first.Data(), man->GetCurrentPeriod().Data(), pass));
   *  pt3TRD.Draw();
   */
  fCanvas->cd(1);
  SetupPadStyle();
  hLikeP_TRD_3tls_electron->Draw("colz");
  fCanvas->cd(2);
  SetupPadStyle();
  hLikeP_TRD_3tls_pion->Draw("colz");
  fCanvas->cd(3);
  SetupPadStyle();
  hLikeP_TRD_4tls_electron->Draw("colz");
  fCanvas->cd(4);
  SetupPadStyle();
  hLikeP_TRD_4tls_pion->Draw("colz");
  fCanvas->cd(5);
  SetupPadStyle();
  hLikeP_TRD_5tls_electron->Draw("colz");
  fCanvas->cd(6);
  SetupPadStyle();
  hLikeP_TRD_5tls_pion->Draw("colz");

  fCanvas->Update();
  fCanvas->Clear();
//.........这里部分代码省略.........

void plot_efficiencies( TFile* file, Int_t type = 2, TDirectory* BinDir=0)
{
   // input:   - Input file (result from TMVA),
   //          - type = 1 --> plot efficiency(B) versus eff(S)
   //                 = 2 --> plot rejection (B) versus efficiency (S)
   //                 = 3 --> plot 1/eff(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 if (type == 3) {
      y1 = 0;
      y2 = -1; // will be set to the max found in the histograms

   } 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";
   if (type == 3) ytit = "1/(Background eff.)";
   TString ftit = ytit + " versus " + xtit;

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


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

   TList xhists;
   TList xmethods;
   UInt_t xnm = TMVAGlob::GetListOfMethods( xmethods );
   TIter xnext(&xmethods);
   // loop over all methods
   TKey *xkey;
   while ((xkey = (TKey*)xnext())) {
      TDirectory * mDir = (TDirectory*)xkey->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() );
         TKey *hkey2;
         while ((hkey2 = TMVAGlob::NextKey(nextKey,"TH1"))) {
            TH1 *h = (TH1*)hkey2->ReadObj();
            TString hname = h->GetName();
            if (hname.Contains( hNameRef ) && hname.BeginsWith( "MVA_" )) {
               if (type==3 && h->GetMaximum() > y2) y2 = h->GetMaximum();
            }
         }
      }
   }


   // 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 );
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	legDta->AddEntry(gr,f->label,"l");
      }
    }
    f++;
  }
    

  mgrDta->Draw("p");
  //  legDta->Draw();   don't draw the data legend
  
  TMultiGraph* mgrThry = new TMultiGraph("Theory","G_{E}^{n}");
  TLegend *legThry = new TLegend(.54,.6,.875,.9,"","brNDC");

  wgr = mgrThry;
  wlg = legThry;

  // the theory
  wlg->SetBorderSize(0); // turn off border
  wlg->SetFillStyle(0);
  
  f = theoryfiles1;
  gr=0;
  while ( f && f->filename ) {
    gr=OneGraph(f);
    if (gr) {
      TGraphAsymmErrors *egr = dynamic_cast<TGraphAsymmErrors*>(gr);
      if (egr && egr->GetN()>1 && egr->GetEYhigh() && egr->GetEYhigh()[1]>0) {
	gr = toerror_band(egr);
	gr->SetFillStyle(3000+f->style);
      }
      if (f->lnpt) {
	wgr->Add(gr,f->lnpt);
	wlg->AddEntry(gr,f->label,f->lnpt);
      }
      else if (gr->GetMarkerStyle()>=20) {
	wgr->Add(gr,"p");
	wlg->AddEntry(gr,f->label,"p");
      }	
      else {
	wgr->Add(gr,"l");
	wlg->AddEntry(gr,f->label,"l");
      }
    }
    f++;
  }

  genf->Draw("same");
  mgrThry->Draw("c");
  legThry->AddEntry(genf,"F_{2}/F_{1} #propto ln^{2}(Q^{2}/#Lambda^{2})/Q^{2}","l");
  legThry->Draw();

  // draw a line at 1
  cn->Modified();

  cn->Update();
  cn->SaveAs(Form("%s.eps",psfile));
  cn->SaveAs(Form("%s.root",psfile));
  gSystem->Exec(Form("./replace_symbols.pl %s.eps",psfile));
  // now an overlay, hopefully matching dimensions

  // remove everything but the graph
  cn->Update();
  TList *clist = cn->GetListOfPrimitives();
  TFrame* frame = cn->GetFrame();
  for (int i=0; i<clist->GetSize(); ) {
    if (clist->At(i) != frame) {
      clist->RemoveAt(i);
    } else i++;
  }
  // draw markers in the corners
  TMarker *mkr = new TMarker(frame->GetX1(),frame->GetY1(),2);
  mkr->Draw();
  mkr = new TMarker(frame->GetX2(),frame->GetY1(),2);
  mkr->Draw();
  mkr = new TMarker(frame->GetX1(),frame->GetY2(),2);
  mkr->Draw();
  mkr = new TMarker(frame->GetX2(),frame->GetY2(),2);
  mkr->Draw();
  frame->SetLineColor(10);
  cn->Update();

  datafile_t miller = { "figure_input/Miller/lattice.GEn.rtf","Miller",
			"[0]","[1]","[1]-[3]","[2]-[1]",0,0,1,3,"F" };

  gr = OneGraph(&miller);
  TGraphAsymmErrors* egr = dynamic_cast<TGraphAsymmErrors*>(gr);
  if (egr && egr->GetEYhigh() && egr->GetEYhigh()[egr->GetN()/2]>0) {
    gr = toerror_band(egr);
    gr->SetLineStyle(1);
    gr->SetFillColor(gr->GetLineColor());
    gr->SetFillStyle(3000+miller.style);
  }
  
  gr->Draw("F");

  cn->Update();
  cn->SaveAs("gen_Miller_Overlay.eps");
  cn->SaveAs("gen_Miller_Overlay.root");
  
}

///
/// Process a directory recursively.
///
void html_a_directory(TDirectory *f, TString path, TEnv *params)
{
  TCanvas *c_h1 = 0;
  if (c_h1 == 0) {
    int x = params->GetValue("H1.XSize", 150);
    int y = params->GetValue("H1.YSize", 100);
    c_h1 = new TCanvas ("c_h1", "1d plots", x, y);
  }

  ///
  /// Check how to make gif plots
  ///

  char command[512];
  sprintf(command, "which pstoimg &> /dev/null");
  bool UsePstoimg = ! system(command);

  ///
  /// Generate the output directory
  ///

  gSystem->MakeDirectory (path);

  ///
  /// Get the html started
  ///

  ofstream html (path + "/index.html");
  html << "<html><head><title>" << f->GetName() << "</title></head>" << endl;
  html << "<body>" << endl;
  html << "<h1>" << f->GetName() << "</h1>" << endl;
  cout << "Processing directory " << f->GetName() << endl;

  ///
  /// Now loop over all the keys in the directory
  ///

  f->cd();
  TList *objlist = f->GetListOfKeys();
  objlist->Sort(); // order alphabetically, instead of order in which they were written
  TIterator *itr = objlist->MakeIterator();
  TKey *key;
  while ((key = static_cast<TKey*>(itr->Next())) != 0) {
    TObject *obj = key->ReadObj();

    if (obj->IsA()->InheritsFrom("TDirectory")) {
      TDirectory *d = static_cast<TDirectory*>(obj);
      html << "<br><a href=\"" << d->GetName() << "/\">" << d->GetName() << "</a>" << endl;
      html_a_directory(d, path + "/" + d->GetName(), params);
    }

    else if (obj->IsA()->InheritsFrom("TObjString")) {
      TObjString *s = static_cast<TObjString*>(obj);
      html << "<p><h2>" << key->GetName() << "</h2>" << endl;
      //html << "<blockquote><pre>" << static_cast<char*>(s->GetString())
	  // << "</pre></blockquote></p>"
	  // << endl;
      html << "<blockquote><pre>" << (s->GetString()).Data() << "</pre></blockquote></p>"<< endl;
      

    }

    //    else if (obj->IsA()->InheritsFrom("CutFlowTable")) {
    //      CutFlowTable *c = static_cast<CutFlowTable*> (obj);
    //      
    //      html << "<p><h2>" << key->GetName() << "</h2>" << endl;
    //
    //      CFTPrinterHTML txt (html);
    //      f->cd();
    //      c->PrintGlobal (txt, "All Events", "");
    //      html << "</p>" << endl;
    //    }

    else if (obj->IsA()->InheritsFrom("TCanvas")) {
      TCanvas *cnv = static_cast<TCanvas*>(obj);

      cnv->Draw();
      cnv->SaveAs(path + "/" + key->GetName() + ".eps");
      if (UsePstoimg) {
	sprintf(command, "pstoimg -type=gif %s &> /dev/null",(path + "/" + key->GetName() + ".eps").Data());
	if (system(command) != 0) {
	  cout<<"Could not convert to gif: "<<path + "/" + key->GetName() + ".eps"<<endl;
	  abort();
	}
      } else cnv->SaveAs(path + "/" + key->GetName() + ".gif");

      cnv->Close();

      html << "<p><a href=\"" << key->GetName() << ".eps\">";
      html << "<img src=\"" << key->GetName() << ".gif\">";
      html << "</a> <br> " << key->GetName() << ".gif </p>" << endl;								   
    }

    else if (obj->IsA()->InheritsFrom("TH1") && !(obj->IsA()->InheritsFrom("TH2"))) {
      TH1 *h = static_cast<TH1*> (obj);

      c_h1->cd();
//.........这里部分代码省略.........

int main(int argc, char *argv[]) 
{
  if (argc != 3 && argc != 4 && argc != 7) {
    cout << "Usage: quickDraw path drawstring [selection [nbins xmin xmax]]" << endl;
    cout << "       will create a file histogram.pdf from the given files" << endl 
	 << endl;
    cout << "path           points to a directory containing ROOT files" << endl;
    cout << "drawstring     is a string usable in TTree::Draw()" << endl;
    cout << "selection      e.g. (muo_n>0)&&(jet_n>0) or 1. if no selection" << endl;
    cout << "nbins          number of bins in histogram" << endl;
    cout << "xmin           low x edge of histogram" << endl;
    cout << "xmax           high x edge of histogram" << endl;
    cout << endl;
    cout << "Example: quickDraw ~/data \"muo_pt[0]\" \"(jet_n>2)\" 100 0 1000" << endl;
    cout << "          will draw muo_pt of the leading muon (index 0) for all events" << endl;
    cout << "          that contain at least three jets in a histogram with 100 bins," << endl;
    cout << "          ranging from 0 to 1000 GeV" << endl;
    cout << endl;
    cout << "Example: quickDraw ~/data vtx_n global_weight" << endl;
    cout << "         will draw vtx_n for all events in an automatically binned histogram " << endl;
    cout << "         The histogram weight will be taken from the leaf global_weight for each event" << endl;
    exit (1);
  }
  // Loop over all files in directory <path> and fill
  // root file names into list <filelist>
  Text_t   lpath[strlen(argv[1])+8];
  sprintf(lpath, "%s", argv[1]);
  Text_t * basepath  = gSystem->ExpandPathName(lpath); // delete it later...
  void   * dirhandle = gSystem->OpenDirectory(basepath);
  const Text_t * basename;
  if (gLogLevel)
    cout << "Reading file names from directory " << basepath << endl;
  TList  * filelist  = new TList;
  while ((basename = gSystem->GetDirEntry(dirhandle))) {
    // Skip non-ROOT files
    if (!strstr(basename, ".root")) 
      continue;
    Text_t * fullname = new Text_t[strlen(basepath)+strlen(basename)+2];
    strcpy(fullname, basepath);
    strcat(fullname, "/");
    strcat(fullname, basename);
    filelist->Add(new TObjString(fullname));
    delete fullname;
  }

  // check if we have read some files
  if (filelist->GetSize() < 1) {
    cerr << "ERR: No *.root files available in " << basepath << " Exiting. " << endl;
    exit(EXIT_FAILURE);
  }

  // Add all root files to a TChain
  if (gLogLevel)
    cout << "Adding files to TChain" << endl;

  TChain * chain = new TChain("ACSkimAnalysis/allData");
  TIter next(filelist);
  while (TObjString * obj = (TObjString *) next()) {
    chain->Add(obj->GetString());
    if (gLogLevel > 3)
      cout << "FILE: " << obj->GetString() << " added to TChain" << endl;
  }

  // Delete all filelist elements and the list itself
  // (They are not deleted automatically when the TList is deleted!!!)
  filelist->Delete();
  delete filelist;
  delete basepath;

  long int events_from_tree = chain->GetEntries();
  cout << "Found a total of " << events_from_tree << " events in chain" << endl;
  
  gROOT->SetBatch();
  TCanvas * c1 = new TCanvas();
  TFile * f  = new TFile("histo.root", "RECREATE");
  Int_t nbins = 200;
  double xmin = 0;
  double xmax = 1;
  if (argc == 7) {
    nbins = atoi(argv[4]);
    xmin = atof(argv[5]);
    xmax = atof(argv[6]);
  }
  TH1D * h1 = new TH1D("h1", argv[2], nbins, xmin, xmax);
  if (argc < 5) {
    h1->SetBit(TH1::kCanRebin);
  }

  const char * selection = "1.";
  if (argc >= 4)
    selection = argv[3];
  
  chain->Draw(Form("%s>>h1", argv[2]), selection);
  h1->Draw();
  c1->Print("histo.pdf");
  f->Write();
  f->Close();
  delete f;
  return 0;
}

void checkRsnPIDqa(TString filename, TString foldername, Bool_t savePng,
		   TString plotTPCpi, TString plotTPCka, TString plotTPCpro,
		   TString plotTTOFpi, TString plotTOFka, TString plotTOFpro)
{
  //Open input file 
  TFile * fin = TFile::Open(filename.Data());
  if (!fin) return 0x0;
  
  //Access output of specific wagon
  TList * list = (TList*) fin->Get(foldername.Data());
  if (!list) return 0x0;

  //Set range for fit
  Float_t RangeFitMomMin = 0.1; //range in momentum where to check the mean and pull
  Float_t RangeFitMomMax = 2.0;
  Int_t xbinFitMin = 0;
  Int_t xbinFitMax = -1;
  Float_t RangeFitNsigmaPIDmin = -2.0; //range in Nsigma where the fit is to be performed
  Float_t RangeFitNsigmaPIDmax = 2.0;

  //Set range for visualisation
  Float_t RangeShowTPC[2] = {0.1, 2.0}; 
  Float_t RangeShowTOF[2] = {0.25, 2.0};
    
  //--------------------------
  // TPC PID Nsigma
  // fit with simple gaussian
  //--------------------------
  //Gaussian function
  TF1 *fGaus = new TF1("f","gaus", -7.0, 7.0);

  //--- pions
  TH2F * hTPCsigmaPi = (TH2F*)list->FindObject(plotTPCpi.Data());
  hTPCsigmaPi->RebinX(2);
  hTPCsigmaPi->SetTitle("TPC Pions");
  MakeUpHisto(hTPCsigmaPi,"p_{TPC} (GeV/c)", "N#sigma_{TPC}", 1, kBlack, 2);
  hTPCsigmaPi->GetYaxis()->SetRangeUser(-5.1,5.1);
  hTPCsigmaPi->GetXaxis()->SetRangeUser(RangeShowTPC[0], RangeShowTPC[1]);
  xbinFitMin = hTPCsigmaPi->GetXaxis()->FindBin(RangeFitMomMin);
  xbinFitMax = hTPCsigmaPi->GetXaxis()->FindBin(RangeFitMomMax);
  hTPCsigmaPi->FitSlicesY(fGaus, xbinFitMin, xbinFitMax );
  TH1D * hTPCsigmaPi_mean = ((TH1D*)gDirectory->FindObject(Form("%s_1", plotTPCpi.Data())))->Clone("hNsigmaTPCpi_mean");
  TH1D * hTPCsigmaPi_pull = ((TH1D*)gDirectory->FindObject(Form("%s_2", plotTPCpi.Data())))->Clone("hNsigmaTPCpi_pull");
  MakeUpHisto(hTPCsigmaPi_mean, "", "", 1, kBlack, 2);
  MakeUpHisto(hTPCsigmaPi_pull, "", "", 1, kRed+2, 2);

  //--- kaons
  TH2F * hTPCsigmaKa = (TH2F*)list->FindObject(plotTPCka.Data());
  hTPCsigmaKa->RebinX(2);
  hTPCsigmaKa->SetTitle("TPC Kaons");
  hTPCsigmaKa->GetYaxis()->SetRangeUser(-5.1,5.1);
  hTPCsigmaKa->GetXaxis()->SetRangeUser(RangeShowTPC[0], RangeShowTPC[1]);
  hTPCsigmaKa->FitSlicesY(fGaus, xbinFitMin, xbinFitMax );
  MakeUpHisto(hTPCsigmaKa,"p_{TPC} (GeV/c)", "N#sigma_{TPC}", 1, kBlack, 2);
  TH1D * hTPCsigmaKa_mean = ((TH1D*)gDirectory->FindObject(Form("%s_1", plotTPCka.Data())))->Clone("hNsigmaTPCka_mean");
  TH1D * hTPCsigmaKa_pull = ((TH1D*)gDirectory->FindObject(Form("%s_2", plotTPCka.Data())))->Clone("hNsigmaTPCka_pull");
  MakeUpHisto(hTPCsigmaKa_mean, "", "", 1, kBlack, 2);
  MakeUpHisto(hTPCsigmaKa_pull, "", "", 1, kRed+2, 2);

  //--- protons
  TH2F * hTPCsigmaPro = (TH2F*)list->FindObject(plotTPCpro.Data());
  hTPCsigmaPro->RebinX(2);
  hTPCsigmaPro->SetTitle("TPC Protons");
  MakeUpHisto(hTPCsigmaPro,"p_{TPC} (GeV/c)", "N#sigma_{TPC}", 1, kBlack, 2);
  hTPCsigmaPro->GetYaxis()->SetRangeUser(-5.1,5.1);
  hTPCsigmaPro->GetXaxis()->SetRangeUser(RangeShowTPC[0], RangeShowTPC[1]);
  hTPCsigmaPro->FitSlicesY(fGaus, xbinFitMin, xbinFitMax );
  TH1D * hTPCsigmaPro_mean = ((TH1D*)gDirectory->FindObject(Form("%s_1", plotTPCpro.Data())))->Clone("hNsigmaTPCpro_mean");
  TH1D * hTPCsigmaPro_pull = ((TH1D*)gDirectory->FindObject(Form("%s_2", plotTPCpro.Data())))->Clone("hNsigmaTPCpro_pull");
  MakeUpHisto(hTPCsigmaPro_mean, "", "", 1, kBlack, 2);
  MakeUpHisto(hTPCsigmaPro_pull, "", "", 1, kRed+2, 2);

   //--- plot TPC
  TLine *l11=new TLine(RangeShowTPC[0],0.,RangeShowTPC[1],0.); l11->SetLineWidth(1); l11->SetLineStyle(7);
  TLine *l12=new TLine(RangeShowTPC[0],1.,RangeShowTPC[1],1.); l12->SetLineWidth(1); l12->SetLineStyle(7);

  gStyle->SetOptStat(0);
  TCanvas *cPidPerformance4 = new TCanvas("cPIDperformance4","TPC PID",1200,500);
  cPidPerformance4->Divide(3,1);
  cPidPerformance4->cd(1);
  gPad->SetLogz(); gPad->SetLogx(); gPad->SetGridx(); gPad->SetGridy();
  hTPCsigmaPi->DrawCopy("colz");
  hTPCsigmaPi_mean->DrawCopy("same");
  hTPCsigmaPi_pull->DrawCopy("same");
  l11->Draw("same"); l12->Draw("same");

  cPidPerformance4->cd(2);
  gPad->SetLogz(); gPad->SetLogx(); gPad->SetGridx(); gPad->SetGridy();
  hTPCsigmaKa->DrawCopy("colz");
  hTPCsigmaKa_mean->DrawCopy("same");
  hTPCsigmaKa_pull->DrawCopy("same");
  l11->Draw("same"); l12->Draw("same");

  cPidPerformance4->cd(3);
  gPad->SetLogz(); gPad->SetLogx(); gPad->SetGridx(); gPad->SetGridy();
  hTPCsigmaPro->DrawCopy("colz");
  hTPCsigmaPro_mean->DrawCopy("same");
  hTPCsigmaPro_pull->DrawCopy("same");
  l11->Draw("same"); l12->Draw("same");

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

void browseStacks( bool makePictures=false, bool wait=true , bool addHistName = false, Double_t maxYScaleF = 1., 
                  bool logScale = false, bool setMinZero = true) {


  gStyle->SetOptTitle(0);

  bool keep2D=false;

  //fix the hNJet histos
  TList *list = gDirectory->GetList();
  TIterator *iter = list->MakeIterator();
  TObject *obj = 0;
  while(obj = iter->Next()) {
  
    if(TString(obj->GetName()).Contains("hnJet") && obj->InheritsFrom(TH1::Class())) {
      
      int nbins = ((TH1F*)obj)->GetNbinsX();
      float overflow = ((TH1F*)obj)->GetBinContent(nbins+1);
      float lastbinval = ((TH1F*)obj)->GetBinContent(nbins);
      ((TH1F*)obj)->SetBinContent(nbins, overflow+lastbinval);
      ((TH1F*)obj)->GetXaxis()->SetBinLabel(nbins, "#geq4");
    }
  }
  
    
    
    
  // Find out what the names of the existing histograms are
  // The histogram names are XX_YY_ZZ, where XX is the sample,
  // eg, "tt", YY is the actual name, ZZ is the final state, eg, "ee"
  TObjArray* myNames = getMyHistosNames("ttdil","ee",keep2D);
    

  // Now loop over histograms, and make stacks
  TCanvas *c = new TCanvas();
  c->Divide(2,2);
  char* suffix[4];
  suffix[0] = "ee";
  suffix[1] = "mm";
  suffix[2] = "em";
  suffix[3] = "all";
  if (makePictures) c->Print("out/stacks.ps[");
  for (int i=0; i<myNames->GetEntries(); i++) {
     
    for (int sample=0; sample<4; sample++) {
       
       
      hist::stack(Form("st_%s_%s",myNames->At(i)->GetName(),suffix[sample]),
		  Form("%s_%s$",myNames->At(i)->GetName(), suffix[sample]));
      THStack* thisStack = (THStack*) gROOT->FindObjectAny(
							   Form("st_%s_%s", myNames->At(i)->GetName(), suffix[sample]));
       
      thisStack->SetMaximum(thisStack->GetMaximum()*maxYScaleF);
      if(TString(myNames->At(i)->GetName()).Contains("hnJet")) {
	TList* histolist = thisStack->GetHists();
	int hatchcount = 0;
	// 	for(int j = 0; j<histolist->GetSize();j++) {
	// 	  if(TString(histolist->At(j)->GetName()).Contains("tt") ||
	// 	     TString(histolist->At(j)->GetName()).Contains("tautau") ||
	// 	     TString(histolist->At(j)->GetName()).Contains("ww") ) continue;
	// 	  hatch(histolist->At(j)->GetName(), FavoriteHatches[hatchcount]);
	// 	  hatchcount++;
	// 	}
      }
	 
	 
      TLegend* thisLeg = hist::legend(thisStack, "lpf", 0, 0, 0.75, 0.65, 0.99, 0.99);
      c->cd(sample+1);
      if (logScale) gPad->SetLogy(); else gPad->SetLogy(0);
      double stackMax = ((TH1*)thisStack->GetHists()->At(0))->GetMaximum();
      double stackMin = ((TH1*)thisStack->GetHists()->At(0))->GetMinimum();
      thisStack->SetMinimum(stackMin);
      if (setMinZero) thisStack->SetMinimum(0);
      if (logScale && stackMin <=0) thisStack->SetMinimum(1e-2*stackMax);
      if (logScale && stackMax == 0) thisStack->SetMinimum(1e-12); 
      thisStack->Draw("hist");
      string xtitle( ((TH1*)gROOT->FindObjectAny(Form("ttdil_%s_%s", myNames->At(i)->GetName(), suffix[sample])))->GetXaxis()->GetTitle());
      string ytitle( ((TH1*)gROOT->FindObjectAny(Form("ttdil_%s_%s", myNames->At(i)->GetName(), suffix[sample])))->GetYaxis()->GetTitle());
      thisStack->GetXaxis()->SetTitle(xtitle.c_str());
      thisStack->GetYaxis()->SetTitle(ytitle.c_str());
      TString hname = thisStack->GetName();
      if(hname.Contains("hnJet")) {
	thisStack->GetXaxis()->SetLabelSize(0.075);
	thisStack->GetYaxis()->SetLabelSize(0.05);
	thisStack->GetXaxis()->SetTitle("N_{jets}");
      }
      thisLeg->Draw();
	
      TPaveText *pt1 = new TPaveText(0.1, 0.95, 0.4, 0.999, "brNDC");
      pt1->SetName("pt1name");
      pt1->SetBorderSize(0);
      pt1->SetFillStyle(0);
	
      TText *blah;
      if (addHistName) blah = pt1->AddText(hname);
      else blah = pt1->AddText("CMS Preliminary");
      blah->SetTextSize(0.05);
      pt1->Draw();
      c->Modified();
	
//.........这里部分代码省略.........

void new_TrackValHistoPublisher(const char* newFile="NEW_FILE",const char* refFile="REF_FILE") {

  cout << ">> Starting new_TrackValHistoPublisher(" 
       << newFile << "," << refFile << ")..." << endl;

  //====  To be replaced from python ====================
  
  const char* dataType = "DATATYPE";
  const char* refLabel("REF_LABEL, REF_RELEASE REFSELECTION");
  const char* newLabel("NEW_LABEL, NEW_RELEASE NEWSELECTION");


  // ==== Initial settings and loads
  gROOT ->SetBatch();
  gErrorIgnoreLevel = kWarning; // Get rid of the info messages
  SetGlobalStyle();

  // ==== Some cleaning... is this needed?  
  delete gROOT->GetListOfFiles()->FindObject(refFile);
  delete gROOT->GetListOfFiles()->FindObject(newFile); 
  
  // ==== Opening files, moving to the right branch and getting the list of sub-branches
  cout << ">> Opening files, moving to the right branch and getting the list of sub-branches..." << endl;

  cout << ">> Finding new DQM file ..." << endl;
  TFile * sfile = new TFile(newFile);
  TList* sl = GetListOfBranches(dataType, sfile);

  if (!sl) {
    cout << "ERROR: Could not find keys!!!" << endl;
    cerr << "ERROR: Could not find keys!!!" << endl;
    return;
  }
  TDirectory*  sdir  = gDirectory;

  if (DEBUG) {
    for (unsigned int i = 0; i < sl->GetEntries(); i++)
      cout << "   + " << sl->At(i)->GetName() << endl;
  }
  
  cout << ">> Finding reference DQM file ..." << endl;
  TFile * rfile = new TFile(refFile);
  TList* rl = GetListOfBranches(dataType, rfile);

  if (!rl) {
    cout << "ERROR: Could not find keys!!!" << endl;
    cerr << "ERROR: Could not find keys!!!" << endl;
    return;
  }
  TDirectory* rdir  = gDirectory;

  if (DEBUG) {
    for (unsigned int i = 0; i < rl->GetEntries(); i++)
      cout << "   + " << rl->At(i)->GetName() << endl;
  }

  //==== Iterate now over histograms and collections
  cout << ">> Iterating over histograms and collections..." << endl;

  bool logy[6]     = {false,  false,  false,  false,   false,  false  };
  bool logx[6]     = {false,  false,  false,  false,   false,  false  };
  bool doKolmo[6]  = {true,   true,   true,   true,    true,   true };
  Double_t norm[6] =  {-1.,-1.,-1.,-1.,-1.,-1.};  // initial default: do not normalize
  Double_t minx[6] = {0, 0, 0, 0, 0, 0};
  Double_t maxx[6] = {0, 0, 0, 0, 0, 0};
  Double_t miny[6] = {0, 0, 0, 0, 0, 0};
  Double_t maxy[6] = {0, 0, 0, 0, 0, 0};
  const char* drawopt[6] = {"", "", "", "", "", ""};
  TString plots[6]       = {"", "", "", "", "", ""};
  TString titles[6]      = {"", "", "", "", "", ""};

  TString rcollname;
  TString scollname;
  TIter iter_r( rl );
  TIter iter_s( sl );
  TString newDirBase("NEW_RELEASE/NEWSELECTION/NEW_LABEL/");
  TKey* rKey = 0;

  // before CMSSW_10_1_0_pre1 a few collection names were different
  bool NEWcollNames = false;
  TString Ref_CMSSW_Release("REF_RELEASE");
  if (Ref_CMSSW_Release.Contains("CMSSW_9") || Ref_CMSSW_Release.Contains("CMSSW_10_0")) NEWcollNames=true;

  while ( (rKey = (TKey*)iter_r()) ) {
    TString myName = rKey->GetName();
    rcollname = myName;
    if (DEBUG) {
      cout << " Checking collection: " << myName << endl;
      cerr << " Checking collection: " << myName << endl;
    }

    TString myName2 = myName;
    if (NEWcollNames) {
      if (myName=="NEWprobeTrks") myName2="probeTrks";
      else if (myName=="NEWprobeTrks_TkAsso") myName2="probeTrks_TkAsso";
      else if (myName=="NEWseedsOfSTAmuons") myName2="seedsOfSTAmuons";
      els