//______________________________________________________________________________
void AODmerge()
{
  // Merging method. No staging and no terminate phase.
  TStopwatch timer;
  timer.Start();
  TString outputDir = "wn.xml";
  TString outputFiles = VAR_AOD_MERGE_FILES;
  TString mergeExcludes = "";
  TObjArray *list = outputFiles.Tokenize(",");
  TIter *iter = new TIter(list);
  TObjString *str;
  TString outputFile;
  Bool_t merged = kTRUE;
  while((str=(TObjString*)iter->Next())) {
    outputFile = str->GetString();
    // Skip already merged outputs
    if (!gSystem->AccessPathName(outputFile)) {
      printf("Output file <%s> found. Not merging again.",outputFile.Data());
      continue;
    }
    if (mergeExcludes.Contains(outputFile.Data())) continue;
    merged = AliAnalysisAlien::MergeOutput(outputFile, outputDir, 10, 0);
    if (!merged) {
      printf("ERROR: Cannot merge %s\n", outputFile.Data());
      return;
    }
  }
  // all outputs merged, validate
  ofstream out;
  out.open("outputs_valid_merge", ios::out);
  out.close();
  timer.Print();
}

 size_t fill_polygon_unique(TIter it, TIter end) {
     size_t num_points = 0;
     osmium::Location last_location;
     for (; it != end; ++it) {
         if (last_location != it->location()) {
             last_location = it->location();
             m_impl.polygon_add_location(m_projection(last_location));
             ++num_points;
         }
     }
     return num_points;
 }

void InitChannels(TTree *tree = (TTree*)gROOT->FindObjectAny("FragmentTree")) {
   if(!tree)
      return;
   TList *list = tree->GetUserInfo();
   TIter *iter = new TIter(list);
   TChannel *chan = 0;
   while(chan = (TChannel*)iter->Next()) {
      TChannel *cur_chan = TChannel::GetChannel(chan->GetAddress());
      TChannel::CopyChannel(cur_chan,chan);
   }

}

 int fill_linestring_unique(TIter it, TIter end) {
     int num_points = 0;
     osmium::Location last_location;
     for (; it != end; ++it) {
         if (last_location != it->location()) {
             last_location = it->location();
             m_impl.linestring_add_location(m_projection(last_location));
             ++num_points;
         }
     }
     return num_points;
 }

void QAmerge(const char *suffix, const char *dir, Int_t stage)
{
// Merging method
  TStopwatch timer;
  timer.Start();
  TString outputDir = dir;
  TString outputFiles = Form("QAresults%s.root,EventStat_temp%s.root,RecoQAresults%s.root",suffix,suffix,suffix);
  TString mergeExcludes = "";
  TObjArray *list = outputFiles.Tokenize(",");
  TIter *iter = new TIter(list);
  TObjString *str;
  TString outputFile;
  Bool_t merged = kTRUE;
  while((str=(TObjString*)iter->Next())) {
    outputFile = str->GetString();
    // Skip already merged outputs
    if (!gSystem->AccessPathName(outputFile)) {
       printf("Output file <%s> found. Not merging again.",outputFile.Data());
       continue;
    }
    if (mergeExcludes.Contains(outputFile.Data())) continue;
    merged = IlcAnalysisAlien::MergeOutput(outputFile, outputDir, 10, stage);
    if (!merged && !outputFile.Contains("RecoQAresults")) {
       printf("ERROR: Cannot merge %s\n", outputFile.Data());
       return;
    }
  }
  // read the analysis manager from file
  if (!outputDir.Contains("Stage")) {
    ofstream out;
    out.open("outputs_valid", ios::out);
    out.close();    
    return;
  }
  IlcAnalysisManager *mgr = IlcAnalysisManager::GetAnalysisManager();
  mgr->SetRunFromPath(mgr->GetRunFromAlienPath(dir));
  mgr->SetSkipTerminate(kFALSE);
  if (!mgr->InitAnalysis()) return;
  mgr->PrintStatus();
  IlcLog::SetGlobalLogLevel(IlcLog::kError);
  TTree *tree = NULL;
  mgr->StartAnalysis("gridterminate", tree);
  if (strlen(suffix)) {
     if (gSystem->Exec(Form("mv trending.root trending%s.root", suffix)))
        ::Error("QAmerge", "File trending.root was not produced");
     if (gSystem->Exec(Form("mv event_stat.root event_stat%s.root", suffix)))
        ::Error("QAmerge", "File trending.root was not produced");
  }   
  ofstream out;
  out.open("outputs_valid", ios::out);
  out.close();
  timer.Print();
}

void QAmerge(const char *dir, Int_t stage)
{
// Merging method
  TStopwatch timer;
  timer.Start();
  TString outputDir = dir;
  TString outputFiles = "QAresults.root,EventStat_temp.root";
  TString mergeExcludes = "";
  TObjArray *list = outputFiles.Tokenize(",");
  TIter *iter = new TIter(list);
  TObjString *str;
  TString outputFile;
  Bool_t merged = kTRUE;
  while((str=(TObjString*)iter->Next())) {
    outputFile = str->GetString();
    // Skip already merged outputs
    if (!gSystem->AccessPathName(outputFile)) {
       printf("Output file <%s> found. Not merging again.",outputFile.Data());
       continue;
    }
    if (mergeExcludes.Contains(outputFile.Data())) continue;
    merged = AliAnalysisAlien::MergeOutput(outputFile, outputDir, 10, stage);
    if (!merged) {
       printf("ERROR: Cannot merge %s\n", outputFile.Data());
       continue;
    }
  }
  TString infolog = "fileinfo.log";
  AliAnalysisAlien::MergeInfo(infolog, dir); 
  // read the analysis manager from file
  if (!outputDir.Contains("Stage")) {
    ofstream out;
    out.open("outputs_valid", ios::out);
    out.close();    
    return;
  }
  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  mgr->SetRunFromPath(mgr->GetRunFromAlienPath(dir));
  mgr->SetSkipTerminate(kFALSE);
  if (!mgr->InitAnalysis()) return;
  mgr->PrintStatus();
  AliLog::SetGlobalLogLevel(AliLog::kError);
  TTree *tree = NULL;
  gROOT->cd();
  mgr->StartAnalysis("gridterminate", tree);
  ofstream out;
  out.open("outputs_valid", ios::out);
  out.close();
  timer.Print();
}

 size_t fill_polygon(TIter it, TIter end) {
     size_t num_points = 0;
     for (; it != end; ++it, ++num_points) {
         m_impl.polygon_add_location(m_projection(it->location()));
     }
     return num_points;
 }

 int fill_linestring(TIter it, TIter end) {
     int num_points = 0;
     for (; it != end; ++it, ++num_points) {
         m_impl.linestring_add_location(m_projection(it->location()));
     }
     return num_points;
 }

//______________________________________________________________________________
void AODmerge()
{
// Merging method. No staging and no terminate phase.
  TStopwatch timer;
  timer.Start();
  TString outputDir = "wn.xml";
  TString outputFiles = "EventStat_temp.root,AODQA.root,AliAOD.root,AliAOD.VertexingHF.root,AliAODGammaConversion.root,FilterEvents_Trees.root,AliAOD.Muons.root";
  TString mergeExcludes = "";
  TObjArray *list = outputFiles.Tokenize(",");
  TIter *iter = new TIter(list);
  TObjString *str;
  TString outputFile;
  Bool_t merged = kTRUE;
  while((str=(TObjString*)iter->Next())) {
    outputFile = str->GetString();
    // Skip already merged outputs
    if (!gSystem->AccessPathName(outputFile)) {
       printf("Output file <%s> found. Not merging again.",outputFile.Data());
       continue;
    }
    if (mergeExcludes.Contains(outputFile.Data())) continue;
    merged = AliAnalysisAlien::MergeOutput(outputFile, outputDir, 10, 0);
    if (!merged) {
       printf("ERROR: Cannot merge %s\n", outputFile.Data());
       continue;
    }
  }
  // all outputs merged, validate
  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  mgr->InitAnalysis();
  mgr->SetGridHandler(new AliAnalysisAlien);
  mgr->StartAnalysis("gridterminate",0);
  ofstream out;
  out.open("outputs_valid", ios::out);
  out.close();
  timer.Print();
}

void draw_network(TDirectory* d)
{
   Bool_t __PRINT_LOGO__ = kTRUE;

   // create canvas
   TStyle* TMVAStyle = gROOT->GetStyle("TMVA"); // the TMVA style
   Int_t canvasColor = TMVAStyle->GetCanvasColor(); // backup
   TMVAStyle->SetCanvasColor( c_DarkBackground );

   static icanvas = -1;
   icanvas++;
   TCanvas* c = new TCanvas( Form( "c%i", icanvas ), Form("Neural Network Layout for: %s", d->GetName()), 
                             100 + (icanvas)*40, 0 + (icanvas+1)*20, 1000, 650  );

   TIter next = d->GetListOfKeys();
   TKey *key;
   TString hName = "weights_hist";
   Int_t numHists = 0;

   // loop over all histograms with hName in name
   while (key = (TKey*)next()) {
      TClass *cl = gROOT->GetClass(key->GetClassName());
      if (!cl->InheritsFrom("TH2F")) continue;    
      TH2F *h = (TH2F*)key->ReadObj();    
      if (TString(h->GetName()).Contains( hName )) {
         numHists++;
      }
   }

   // loop over all histograms with hName in name again
   next.Reset();
   Double_t maxWeight = 0;

   // find max weight
   while (key = (TKey*)next()) {

      //cout << "Title: " << key->GetTitle() << endl;
      TClass *cl = gROOT->GetClass(key->GetClassName());
      if (!cl->InheritsFrom("TH2F")) continue;    

      TH2F* h = (TH2F*)key->ReadObj();    
      if (TString(h->GetName()).Contains( hName )){
 
         Int_t n1 = h->GetNbinsX();
         Int_t n2 = h->GetNbinsY();
         for (Int_t i = 0; i < n1; i++) {
            for (Int_t j = 0; j < n2; j++) {
               Double_t weight = TMath::Abs(h->GetBinContent(i+1, j+1));
               if (maxWeight < weight) maxWeight = weight;
            }
         }
      }
   }

   // draw network
   next.Reset();
   Int_t count = 0;
   while (key = (TKey*)next()) {

      TClass *cl = gROOT->GetClass(key->GetClassName());
      if (!cl->InheritsFrom("TH2F")) continue;    

      TH2F* h = (TH2F*)key->ReadObj();    
      if (TString(h->GetName()).Contains( hName )){
         draw_layer(c, h, count++, numHists+1, maxWeight);
      }
   }

   draw_layer_labels(numHists+1);

   // ============================================================
   if (__PRINT_LOGO__) TMVAGlob::plot_logo();
   // ============================================================  

   c->Update();

   TString fname = "plots/network";
   TMVAGlob::imgconv( c, fname );

   TMVAStyle->SetCanvasColor( canvasColor );
}

void
makeUnfoldingSystematicsComparisonPlots (int whichobservable, int whichjet, int whichlepton, int inclusive)
{
  //Open the file and form the name
  string fileSystematics; 
  string suffix="/gpfs/cms/data/2011/Systematics/postApproval_v58_Journal/";
  //string suffix="/tmp/";

  if (whichlepton==1) suffix=suffix+"ele/";
  if (whichlepton==2) suffix=suffix+"muo/";

  //  setTDRStyle ();
  gStyle->SetErrorX(0);
  gStyle->SetPadGridX(0);
  gStyle->SetPadGridY(0);

  bool absoluteNormalization = true;
  int lepton = 3; //1 -> electron,  2-> muon , 3 -> combined reults!
  bool addLumiUncertainties = true;
  double lumiError = 0.025;
  int use_case = whichobservable;
  int whichjet = whichjet;
  string version = "_v2_32";

  //string s         = "/afs/infn.it/ts/user/marone/html/ZJets/FinalPlotsForAN/v41/SVDBayes/";
  string s="/afs/infn.it/ts/user/marone/html/ZJets/FinalPlotsForAN/v58_Journal/UnfoldingSyst/ele/";
  //  string s         = "/gpfs/cms/users/schizzi/EleMuComparisonPlots/PostUnfolding/";
  string  eleplotpath = "/gpfs/cms/users/schizzi/Systematics/ele/";
  string  muoplotpath = "/gpfs/cms/users/schizzi/Systematics/muo/";

  gStyle->SetOptStat (0);

  TCanvas *plots = new TCanvas ("plots", "EB", 200, 100, 600, 800);

  //DATA:
  //string elepathFile   ="/gpfs/cms/data/2011/Unfolding/testReferenceMu.root";
  //string muopathFile   ="/gpfs/cms/data/2011/Unfolding/testMu.root";
  string elepathFile   ="/gpfs/cms/data/2011/Unfolding/UnfoldingOfficialV58_afterCWR.root";
  string muopathFile   ="/gpfs/cms/data/2011/Unfolding/UnfoldingOfficialV58_BinWidth.root";
  string thirdpathFile   ="/gpfs/cms/data/2011/Unfolding/UnfoldingOfficialV58_SherpaV2.root";

  if (whichlepton==2){
  elepathFile   ="/gpfs/cms/data/2011/Unfolding/UnfoldingOfficialV58_afterCWRMu.root";
  muopathFile   ="/gpfs/cms/data/2011/Unfolding/UnfoldingOfficialV58_BinWidthMu.root";
thirdpathFile   ="/gpfs/cms/data/2011/Unfolding/UnfoldingOfficialV58_SherpaV2Mu.root";
  s="/afs/infn.it/ts/user/marone/html/ZJets/FinalPlotsForAN/v58_Journal/UnfoldingSyst/muo/";
  }
  //string elepathFile   ="/gpfs/cms/data/2011/Unfolding/UnfoldingOfficialV57_3NoSQRT.root";
  //string muopathFile   ="/gpfs/cms/data/2011/Unfolding/UnfoldingOfficialV57_3NoMCToy.root";
    
  TFile *histof = TFile::Open (elepathFile.c_str ());
  histof->cd ("");
  TDirectory *dir = gDirectory;
  TList *mylist = (TList *) dir->GetListOfKeys ();
  TIter iter (mylist);
  TObject *tobj = 0;

  TFile *histof2 = TFile::Open (thirdpathFile.c_str ());
  histof2->cd ("");
  TDirectory *dir2 = gDirectory;
  TList *mylist2 = (TList *) dir2->GetListOfKeys ();
  TIter iter2 (mylist);
  TObject *tobj2 = 0;


  TFile *histofmuo = TFile::Open (muopathFile.c_str ());


  string stringmatch;
  string systPathFile;
  string systPathFileMuo;

  histof->cd ("");

  int i = 0;			// solo di servizio quando debuggo...
  while ((tobj = iter.Next ()))
    {
      string name = tobj->GetName ();

      if (use_case == 1) 
	{ // Jet Multiplicity
	  stringmatch = "JetMultiplicityUnfolded";
	  systPathFile = eleplotpath + "systematicsEff_jetMult" + version + ".txt";
	  systPathFileMuo = muoplotpath + "systematicsEff_jetMult" + version + ".txt";
	  fileSystematics = suffix+"systematicsUnfReweight_jetMult" + version + ".txt"; 
	}

      if (use_case == 2) 
	{ // Jet Pt
	  if (whichjet == 1)
	    {
	      stringmatch = "jReco_leading";
	      systPathFile = eleplotpath + "systematicsEff_jet1Pt" + version + ".txt";
	      systPathFileMuo = muoplotpath + "systematicsEff_jet1Pt" + version + ".txt";
	      fileSystematics = suffix+"systematicsUnfReweight_jet1Pt" + version + ".txt"; 
	    }
	  
	  if (whichjet == 2)
	    {
	      stringmatch = "jReco_subleading";
//.........这里部分代码省略.........

std::string generateCodeFromStreamers(std::string url, std::string treeLocation, std::vector<std::string> &classNames, std::string &errorMessage) {
  TFile *tfile = TFile::Open(url.c_str());
  if (tfile == nullptr  ||  !tfile->IsOpen()) {
    errorMessage = std::string("File not found: ") + url;
    return std::string();
  }

  if (tfile->IsZombie()) {
    errorMessage = std::string("Not a ROOT file: ") + url;
    return std::string();
  }

  TTreeReader reader(treeLocation.c_str(), tfile);
  if (reader.IsZombie()) {
    errorMessage = std::string("Not a TTree: ") + treeLocation.c_str() + std::string(" in file: ") + url;
    return std::string();
  }

  TTree *ttree = reader.GetTree();

  std::set<std::string> includes;
  std::vector<ClassStructure> classes;

  TIter listOfBranches = ttree->GetListOfBranches();
  for (TBranch *tbranch = (TBranch*)listOfBranches.Next();  tbranch != nullptr;  tbranch = (TBranch*)listOfBranches.Next()) {
    TClass *tclass = TClass::GetClass(tbranch->GetClassName());
    if (tclass != nullptr  &&  tbranch->GetListOfBranches()->GetEntries() > 0)
      classesFromBranch(tbranch, tclass, classes, 0, includes);
  }

  for (int i = 0;  i < classes.size();  i++)
    classNames.push_back(classes[i].fullName);

  tfile->Close();

  std::string out;

  for (std::set<std::string>::iterator iter = includes.begin();  iter != includes.end();  ++iter)
    out += *iter + "\n";
  out += "\n";
  
  for (std::vector<ClassStructure>::iterator iter = classes.begin();  iter != classes.end();  ++iter) {
    int i = 0;
    for (;  i < iter->splitName.size() - 1;  i++)
      out += std::string(i * 2, ' ') + "namespace " + iter->splitName[i] + " {\n";

    out += std::string(i * 2, ' ') + "class " + iter->splitName.back() + ";\n";
    i--;
    for (;  i >= 0;  i--)
      out += std::string(i * 2, ' ') + "}\n";
  }
  out += "\n";

  for (std::vector<ClassStructure>::iterator iter = classes.begin();  iter != classes.end();  ++iter) {
    int i = 0;
    for (;  i < iter->splitName.size() - 1;  i++)
      out += std::string(i * 2, ' ') + "namespace " + iter->splitName[i] + " {\n";
    out += iter->cpp(i * 2) + "\n";
    i--;
    for (;  i >= 0;  i--)
      out += std::string(i * 2, ' ') + "}\n";
  }

  for (std::vector<ClassStructure>::iterator iter = classes.begin();  iter != classes.end();  ++iter) {
    int i = 0;
    for (;  i < iter->splitName.size() - 1;  i++)
      out += std::string(i * 2, ' ') + "namespace " + iter->splitName[i] + " {\n";
    out += std::string(i * 2, ' ') + "ClassImp(" + iter->splitName.back() + ")\n";
    i--;
    for (;  i >= 0;  i--)
      out += std::string(i * 2, ' ') + "}\n";
  }
  
  return out;
}

void classesFromBranch(TBranch *tbranch, TClass *tclass, std::vector<ClassStructure> &classes, int prefix, std::set<std::string> &includes) {
  std::string className = tclass->GetName();

  if (className == std::string("TObject"))
    includes.insert("#include \"TObject.h\"");
  else if (className == std::string("TRef"))
    includes.insert("#include \"TRef.h\"");
  else if (className == std::string("TRefArray"))
    includes.insert("#include \"TRefArray.h\"");
  else if (className == std::string("TH1"))
    includes.insert("#include \"TH1.h\"");
  else if (className == std::string("TBits")) {
    includes.insert("#include \"TBits.h\"");
  }
  else {
    bool classSeen = false;
    for (int i = 0;  i < classes.size();  i++)
      if (classes[i].fullName == className)
        classSeen = true;

    if (!classSeen) {
      TVirtualStreamerInfo *tVirtualStreamerInfo = tclass->GetStreamerInfo();
      int version = tVirtualStreamerInfo->GetClassVersion();

      ClassStructure classStructure(tclass, version);

      if (tVirtualStreamerInfo->GetElements()->GetEntries() == 1  &&
          ((TStreamerElement*)(tVirtualStreamerInfo->GetElements()->First()))->IsBase()  &&
          std::string(tVirtualStreamerInfo->GetElements()->First()->GetName()) == std::string("TObjArray")) {

        TVirtualStreamerInfo *substreamer = ((TBranchElement*)tbranch)->GetInfo();

        if (std::string(substreamer->GetName()) == std::string("TClonesArray")) {
          ROOT::Internal::TTreeGeneratorBase ttreeGenerator(tbranch->GetTree(), "");
          TString className = ttreeGenerator.GetContainedClassName((TBranchElement*)tbranch, (TStreamerElement*)(substreamer->GetElements()->First()), true);
          classesFromBranch(tbranch, TClass::GetClass(className), classes, prefix + std::string(tbranch->GetName()).size() + 1, includes);
        }
      }

      TIter elements = tVirtualStreamerInfo->GetElements();
      for (TStreamerElement *tStreamerElement = (TStreamerElement*)elements.Next();  tStreamerElement != nullptr;  tStreamerElement = (TStreamerElement*)elements.Next()) {
        std::string streamerName = tStreamerElement->GetName();

        if (tStreamerElement->IsBase()) {
          TClass *elementClass = tStreamerElement->GetClassPointer();
          std::string type = elementClass->GetName();

          if (type == std::string("TObject")) {
            classStructure.bases.push_back(type);
            includes.insert(std::string("#include \"") + type + std::string(".h\""));
          }
        }

        else {
          TIter listOfBranches = tbranch->GetListOfBranches();
          TBranch *subbranch;
          for (subbranch = (TBranch*)listOfBranches.Next();  subbranch != nullptr;  subbranch = (TBranch*)listOfBranches.Next()) {
            std::string branchName = subbranch->GetName();
            branchName = branchName.substr(prefix, std::string::npos);
            int firstDot = branchName.find('.');
            int firstBracket = branchName.find('[');
            if (firstDot != std::string::npos  &&  firstBracket != std::string::npos)
              branchName = branchName.substr(0, firstDot < firstBracket ? firstDot : firstBracket);
            else if (firstDot != std::string::npos)
              branchName = branchName.substr(0, firstDot);
            else if (firstBracket != std::string::npos)
              branchName = branchName.substr(0, firstBracket);
            if (branchName == streamerName)
              break;
          }
          if (subbranch == nullptr)
            continue;

          std::string branchName = subbranch->GetName();

          std::string variable = tStreamerElement->GetName();
          std::string variableWithArray = subbranch->GetName();
          variableWithArray = variableWithArray.substr(prefix, std::string::npos);
          variableWithArray = variableWithArray.substr(0, variableWithArray.find('.'));

          std::string comment = tStreamerElement->GetTitle();

          std::string type;
          Bool_t pointer = false;
          switch (tStreamerElement->GetType()) {

          case TVirtualStreamerInfo::kOffsetP + TVirtualStreamerInfo::kBool:
            pointer = true;
          case TVirtualStreamerInfo::kOffsetL + TVirtualStreamerInfo::kBool:
          case TVirtualStreamerInfo::kBool:
            type = "Bool_t";
            break;

          case TVirtualStreamerInfo::kOffsetP + TVirtualStreamerInfo::kChar:
            pointer = true;
          case TVirtualStreamerInfo::kOffsetL + TVirtualStreamerInfo::kChar:
          case TVirtualStreamerInfo::kChar:
            type = "Char_t";
            break;

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