C++ TList::Write方法代码示例

void WriteXsection()
{
    TPythia6 *pythia = TPythia6::Instance();
    pythia->Pystat(1);

    Int_t    ntrials  = pythia->GetMSTI(5);
    Double_t xsection = pythia->GetPARI(1);

    cout << "Pythia cross section: " << xsection
         << ", number of trials: " << ntrials << endl;

    TFile *file = TFile::Open("pyxsec_hists.root", "NEW");
    TList *list = new TList();
    TProfile *h1Xsec = new TProfile("h1Xsec", "", 1, 0., 1.);
    list->Add(h1Xsec);
    TH1F  *h1Trials  = new TH1F("h1Trials",   "", 1, 0., 1.);
    h1Trials->Sumw2();
    list->Add(h1Trials);

    h1Xsec->Fill(0.5, xsection);
    h1Trials->Fill(0.5, ntrials);
    file->cd();
    list->Write("cFilterList", TObject::kSingleKey);
    file->Close();

    return;
}

void analysis_t_test1::Terminate()
{
   // The Terminate() function is the last function to be called during
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.
  TFile file_out("analysis_t_test1_output.root","recreate");
  TList *outlist = GetOutputList();
  
  outlist->Write();
  file_out.Write();
  file_out.Close();

}

void neutrootracker::DoAnalysis(const TString filelist, const Int_t ntarget, const TString tag)
{
    gRooTrackerUtils::Ini();
    Init(gRooTrackerUtils::IniIO(filelist,"nRooTracker"));
    printf("NeutTrackerAnalysis ntarget = %d\n", ntarget);
    
    TList *lout = gRooTrackerUtils::GetHistList();
    TTree *tout = gRooTrackerUtils::GetTree();
    
    const Int_t nentries = fChain->GetEntries();
    Int_t nsel = 0;
    Int_t nfail = 0;
    
    for (Int_t jentry = 0; jentry<nentries; jentry++) {
        if(jentry%10000==0) printf("%d/%d nsel %d nfail %d\n", jentry, nentries, nsel, nfail);
        if(!fChain->GetEntry(jentry)){
            printf("DoAnalysis GetEntry ends at jentry %d\n", jentry);
            break;
        }
        
        gRooTrackerUtils::FillCount(lout,0);
        
        const Int_t nfound = gRooTrackerUtils::FindParticles(StdHepN, StdHepPdg, StdHepStatus, StdHepP4, StdHepP4, "NEUT");
        
        gRooTrackerUtils::SetCount(nfound);
        gRooTrackerUtils::SetpTT();
        gRooTrackerUtils::SetDeltapTT();
        gRooTrackerUtils::SetConsSum();
        gRooTrackerUtils::SetAltdpTT();
        gRooTrackerUtils::FillCount(lout,11);
        tout->Fill();
        nsel++;
        if(ntarget>=0 && nsel>=ntarget) break;
    }
    printf("nentries %d nsel %d nfail %d delta %d\n", nentries, nsel, nfail, nentries-nsel-nfail);
    
    if( (nentries-nsel-nfail)!=0 ){
        printf("NeutRooTracker: nsel nfail bad nsel %d nfail %d nentries %d\n", nsel, nfail, nentries);
        //exit(1);
    }
    
    TFile *fout= new TFile(Form("%s_DeltapTT.root",tag.Data()),"recreate");
    tout->Write();
    //for further handling in plot
    lout->Write("lout",TObject::kSingleKey);
    fout->Save();
    fout->Close();
    delete fout;
}

void tables_write()
{
   // first connect to data base
   // "recreate" option delete all your tables !!!!
   TSQLFile* f = new TSQLFile(dbname, "recreate", username, userpass);
   if (f->IsZombie()) { delete f; return; }

   // you can change configuration only until first object
   // is written to TSQLFile
   f->SetUseSuffixes(kFALSE);
   f->SetArrayLimit(1000);
   f->SetUseIndexes(1);
//   f->SetTablesType("ISAM");
//   f->SetUseTransactions(kFALSE);

   // lets first write histogram
   TH1I* h1 = new TH1I("histo1","histo title", 1000, -4., 4.);
   h1->FillRandom("gaus",10000);
   h1->Write("histo");
   h1->SetDirectory(0);

   // here we create list of objects and store them as single key
   // without kSingleKey all TBox objects will appear as separate keys
   TList* arr = new TList;
   for(Int_t n=0;n<10;n++) {
      TBox* b = new TBox(n*10,n*100,n*20,n*200);
      arr->Add(b, Form("option_%d_option",n));
   }
   arr->Write("list",TObject::kSingleKey);

   // clones array is also stored as single key
   TClonesArray clones("TBox",10);
   for(int n=0;n<10;n++)
       new (clones[n]) TBox(n*10,n*100,n*20,n*200);
   clones.Write("clones",TObject::kSingleKey);

   // close connection to database
   delete f;
}

void AdaptFileStructure(TString fileName)
{
  TFile* f = TFile::Open(fileName.Data(), "READ");
  TList* cList = new TList();
  cList->SetName("cList");
  cList->SetOwner(kTRUE);
  for (Int_t i = 0; i < gFile->GetNkeys(); i++) {
    TObject* obj = f->GetListOfKeys()->At(i);
    TString name = obj->GetName();
    TObject* obj2 = f->Get(name.Data());
    cList->Add(obj2); 
  }
  
  TString fileNameOutput = fileName;
  fileNameOutput.ReplaceAll(".root", "_commonStructure.root");
  TFile* fOut = TFile::Open(fileNameOutput.Data(), "RECREATE");
  fOut->mkdir("PWGLF_PPVsMultXCheck_MC");
  fOut->cd("PWGLF_PPVsMultXCheck_MC");
  cList->Write("cList", TObject::kSingleKey);
  fOut->Close();
  f->Close();
}

/**
 * Steering 
 * 
 * @param input Input file 
 */
void
ExtractGSE(const char* input)
{
  if (!gROOT->GetClass("GraphSysErr")) 
    gROOT->LoadMacro("$HOME/GraphSysErr/GraphSysErr.C+g");

  TString base = gSystem->BaseName(input); base.ReplaceAll(".root","");
  TFile*  file = TFile::Open(input, "READ");
  if (!file) return;

  TH1* cent = GetH1(file, "realCent");

  Bool_t first = true;
  TList* stack = new TList;
  for (Int_t i = 1; i <= cent->GetNbinsX(); i++) {
    Double_t c1 = cent->GetXaxis()->GetBinLowEdge(i);
    Double_t c2 = cent->GetXaxis()->GetBinUpEdge(i);
    TObject* g  = MakeGSE(file, c1, c2);
    if (!g) continue;
    stack->Add(g);    
    if (first) g->Draw("quad stat combine axis");
    else       g->Draw("quad stat combine");
    first = false;
  }

  TString obase(base); obase.Prepend("GSE_");
  
  std::ofstream out(Form("%s.input", obase.Data()));
  GraphSysErr::Export(stack, out, "HFC", 2);
  out << "*E" << std::endl;
  out.close();

  TFile* rout = TFile::Open(Form("%s.root", obase.Data()), "RECREATE");
  stack->Write("container", TObject::kSingleKey);
  rout->Write();
  
}

//.........这里部分代码省略.........
   theTree->SetBranchAddress( "deltaRphotonjet", &mydeltaRphotonjet );
   theTree->SetBranchAddress( "deltaRphotonmuon", &mydeltaRphotonmuon );
//   theTree->SetBranchAddress( "ht", &myht );
   theTree->SetBranchAddress( "costopphoton", &mycostopphoton );
   theTree->SetBranchAddress( "jetmultiplicity", &myjetmultiplicity );
//   theTree->SetBranchAddress( "bjetmultiplicity", &mybjetmultiplicity );
   theTree->SetBranchAddress( "deltaphiphotonmet", &mydeltaphiphotonmet );
   theTree->SetBranchAddress( "cvsdiscriminant", &mycvsdiscriminant );
   theTree->SetBranchAddress( "coswphoton", &mycoswphoton );
//   theTree->SetBranchAddress( "leptoncharge", &myleptoncharge );
   theTree->SetBranchAddress( "weight", &myweight);
   theTree->SetBranchAddress( "btagSF", &mybtagSF);
   theTree->SetBranchAddress( "btagSFup", &mybtagSFup);
   theTree->SetBranchAddress( "btagSFdown", &mybtagSFdown);
   theTree->SetBranchAddress( "mistagSFup", &mymistagSFup);
   theTree->SetBranchAddress( "mistagSFdown", &mymistagSFdown);
   theTree->SetBranchAddress( "triggerSF", &mytriggerSF);
   theTree->SetBranchAddress( "triggerSFup", &mytriggerSFup);
   theTree->SetBranchAddress( "triggerSFdown", &mytriggerSFdown);
   theTree->SetBranchAddress( "photonSF", &myphotonSF);
   theTree->SetBranchAddress( "photonSFup", &myphotonSFup);
   theTree->SetBranchAddress( "photonSFdown", &myphotonSFdown);
   theTree->SetBranchAddress( "muonSF", &mymuonSF);
   theTree->SetBranchAddress( "muonSFup", &mymuonSFup);
   theTree->SetBranchAddress( "muonSFdown", &mymuonSFdown);
   theTree->SetBranchAddress( "pileupSF", &mypileupSF);
   theTree->SetBranchAddress( "pileupSFup", &mypileupSFup);
   theTree->SetBranchAddress( "pileupSFdown", &mypileupSFdown);
   theTree->SetBranchAddress( "jetmatchinginfo", &myjetmatchinginfo );

   for (Long64_t ievt=0; ievt<theTree->GetEntries();ievt++) {
   double finalweight;
   theTree->GetEntry(ievt);

finalweight=(*myweight)[0];

//cout<<(*myweight)[0]<<endl;

if((*mymasstop )[0]>mtopdown && (*mymasstop )[0]<mtopup){
Shists[0] ->Fill( (*mydeltaRphotonjet)[0],finalweight*signalscales[phi]);
Shists[1] ->Fill((*mydeltaRphotonjet)[0],finalweight*signalscales[phi]);
Shists[2] ->Fill( (*mydeltaRphotonjet)[0],finalweight*signalscales[phi]);
Shists[3] ->Fill( (*mydeltaRphotonjet)[0],finalweight*signalscales[phi]);
Shists[4] ->Fill( (*mydeltaRphotonjet)[0],finalweight*signalscales[phi]);
}
}

delete myptphoton;
delete myetaphoton;
delete myptmuon;
delete myetamuon;
delete myptjet;
delete myetajet;
delete mymasstop;
//delete mymtw;
delete mydeltaRphotonjet;
delete mydeltaRphotonmuon;
//delete myht;
delete mycostopphoton;
delete mydeltaphiphotonmet;
delete mycvsdiscriminant;
delete myjetmultiplicity;
//delete mybjetmultiplicity;
//delete myleptoncharge;
//delete myplot;
delete mybtagSF;
delete mybtagSFup;
delete mybtagSFdown;
delete mymistagSFup;
delete mytriggerSF;
delete mytriggerSFup;
delete mytriggerSFdown;
delete myphotonSF;
delete myphotonSFup;
delete myphotonSFdown;
delete mypileupSF;
delete mypileupSFup;
delete mypileupSFdown;
delete input;
delete myjetmatchinginfo;
/*Shists[0]->Scale(5/Shists[0]->Integral());
Shists[1]->Scale(10/Shists[1]->Integral());
Shists[2]->Scale(20/Shists[2]->Integral());
Shists[3]->Scale(30/Shists[3]->Integral());
*/
//Shists[4]->Scale(40/Shists[4]->Integral());

//hList->Add(Shists[0]);
//hList->Add(Shists[1]);
//hList->Add(Shists[2]);
//hList->Add(Shists[3]);
hList->Add(Shists[4]);
}


/////////////////////////////////////////////////////////////////////////
TFile *target  = new TFile( "DELTARIPHJET.root","RECREATE" );
  hList->Write();
   target->Close();
}

//.........这里部分代码省略.........
    cerr << "Could not retrieve fileset: " << gFilesetname << endl;
    return;
  }

  vector<string> fileNames = GetDataVector(files,gFfrom,gFto);
  fwlite::ChainEvent event(fileNames);  

  //Initialize the analyzer      
  EbyEFlowAnalyzerFWLite* corr = new EbyEFlowAnalyzerFWLite(event);
  corr->SetCutParameters(gCut);

  TString strtrgid = "Track";
  if(gSystem->Getenv("TRGID")) strtrgid = gSystem->Getenv("TRGID");
  DiHadronCorrelationMultiBaseFWLite::ParticleType gTrgID = corr->GetParticleID(strtrgid);
  TString strassid = "Track";
  if(gSystem->Getenv("ASSID")) strassid = gSystem->Getenv("ASSID");
  DiHadronCorrelationMultiBaseFWLite::ParticleType gAssID = corr->GetParticleID(strassid);
  corr->SetTrgID(gTrgID);
  corr->SetAssID(gAssID);

  if(gCentfilename.Contains("root"))
  {
    cout<<"Running correlation analysis in Heavy Ion!"<<endl;
    TFile* fcentfile = new TFile(gCentfilename.Data());
    if(fcentfile->IsOpen()) corr->SetCentrality(fcentfile,gCenttablename,gNCentBins,gCentRunnum);
    else cout<<"Centrality table cannot be opened!"<<endl;
  }
  else cout<<"Running correlation analysis in pp!"<<endl;

  if(gEffhistname.Contains("root"))
  {
    TFile* feffhist = new TFile(gEffhistname.Data());    
    if(feffhist->IsOpen())
    {
      TH3D* htotaleff = 0;
      if(!gCut.IsHI) htotaleff = (TH3D*)feffhist->Get("rTotalEff3D");
      else
      {
        int histcentbin = -1;
        if((gCut.centmin>=0 && gCut.centmax<=2) || gCut.centmin>=50) histcentbin = 0;
        if(gCut.centmin>=2 && gCut.centmax<=4) histcentbin = 1;
        if(gCut.centmin>=4 && gCut.centmax<=12) histcentbin = 2;
        if(gCut.centmin>=12 && gCut.centmax<=20) histcentbin = 3;
        if(gCut.centmin>=20 && gCut.centmax<=40) histcentbin = 4;
        if(gCut.centmin==-1 && gCut.centmax==-1) histcentbin = 4;

        htotaleff = (TH3D*)feffhist->Get(Form("Tot_%d",histcentbin));
      }
      corr->LoadEffWeight(htotaleff);
      cout<<"Tracking efficiency weighting histogram is loaded!"<<endl;
    }
    else cout<<"Tracking efficiency weighting file cannot be opened!"<<endl;
  }
  else cout<<"No Tracking efficiency weighting histogram is found, or it's running MC!"<<endl;

  if(gTrghistname.Contains("root"))
  {
    TFile* ftrghist = new TFile(gTrghistname.Data());
    if(ftrghist->IsOpen())
    {
      corr->LoadTrgWeight((TH1D*)ftrghist->Get("trgEff"));
      cout<<"Triggering efficiency weighting histogram is loaded!"<<endl;
    }
    else cout<<"Triggering efficiency weighting file cannot be opened!"<<endl;
  }
  else cout<<"No Triggering efficiency weighting histogram is found, or it's running MC!"<<endl;

  if(gPileupdistfunchistname.Contains("root"))
  {
    TFile* fpileupdistfunchist = new TFile(gPileupdistfunchistname.Data());
    if(fpileupdistfunchist->IsOpen())
    {
      corr->LoadPileUpDistFunc((TH1D*)fpileupdistfunchist->Get("distfunc"));
      cout<<"Pileup distfunc histogram is loaded!"<<endl;
    }
    else cout<<"Pileup distfunc file cannot be opened!"<<endl;
  }

  corr->Process();

  // open output file
  TString outputfilename = Form("%s/%s/unmerged/%s_%s_ffrom%d_fto%d_vtxmin%.1f_vtxmax%.1f_nmin%d_nmax%d_etatrg%.1f-%.1f_etaass%.1f-%.1f_centmin%d_centmax%d.root",gOutputDir.Data(),gEventtype.Data(),gEventtype.Data(),gTag.Data(),gFfrom,gFto,gCut.zvtxmin,gCut.zvtxmax,gCut.nmin,gCut.nmax,gCut.etatrgmin,gCut.etatrgmax,gCut.etaassmin,gCut.etaassmax,gCut.centmin,gCut.centmax);
  TFile* outf = new TFile(outputfilename.Data(),"recreate");

  // Save the outputs
  TList* output = corr->GetOutputs();
  if(outf && output) {
    outf->cd();
    output->Write();
    cout << "Output file " <<  outf->GetName()
    << " written: " << endl;
    outf->ls();
    outf->Close();
    if(outf) delete outf;
  }  
  delete output;

  watch->Stop();
  watch->Print();
}

//.........这里部分代码省略.........
    fastjet::ClusterSequenceArea clustSeq(fjInputVac, jetsDef, areaDef);
    std::vector<fastjet::PseudoJet> includJetsPy = clustSeq.inclusive_jets(dJetsPtMin);
//  std::vector<fastjet::PseudoJet> subtedJetsPy = bkgSubtractor(includJetsPy);
    std::vector<fastjet::PseudoJet> selectJetsPy = selectJet(includJetsPy);
//  std::vector<fastjet::PseudoJet> sortedJetsPy = fastjet::sorted_by_pt(selectJetsPy);

    for (int j=0; j<selectJetsPy.size(); j++) {
      SetJetUserInfo(selectJetsPy[j]);
      selectJetsPy[j].set_user_index(j);
    }
//=============================================================================

    bkgsEstimator.set_particles(fjInputHyd);
    double dBkgRhoHd = bkgsEstimator.rho();
    double dBkgRmsHd = bkgsEstimator.sigma();

    fastjet::ClusterSequenceArea clustSeqHd(fjInputHyd, jetsDef, areaDef);
    std::vector<fastjet::PseudoJet> includJetsHd = clustSeqHd.inclusive_jets(dJetsPtMin);
    std::vector<fastjet::PseudoJet> selectJetsHd = selectJet(includJetsHd);

    for (int j=0; j<selectJetsHd.size(); j++) {
      SetJetUserInfo(selectJetsHd[j]);
      selectJetsHd[j].set_user_index(j);
      if (selectJetsHd[j].user_info<UserInfoJet>().IsBkg()) continue;

      for (int i=0; i<selectJetsPy.size(); i++) {
        if (CalcDeltaR(selectJetsHd[j],selectJetsPy[i])>0.8) continue;
        DoTrkMatch(selectJetsHd[j], selectJetsPy[i]);
      }
    }
//=============================================================================

  








    for (int j=0; j<sortedJets.size(); j++) {
      double dJet = sortedJets[j].pt();

      hJet->Fill(dJet, dNorm);
//=============================================================================

      fastjet::Filter trimmer(subjDef, fastjet::SelectorPtFractionMin(0.));
      fastjet::PseudoJet trimmdJet = trimmer(sortedJets[j]);
      std::vector<fastjet::PseudoJet> trimmdSj = trimmdJet.pieces();

      double nIsj = 0.;
      double d1sj = -1.; int k1sj = -1;
      double d2sj = -1.; int k2sj = -1;
      for (int i=0; i<trimmdSj.size(); i++) {
        double dIsj = trimmdSj[i].pt(); if (dIsj<0.001) continue;

        hJetIsj->Fill(dJet, dIsj, dNorm);
        hJetIsz->Fill(dJet, dIsj/dJet, dNorm);

        if (dIsj>d1sj) {
          d2sj = d1sj; k2sj = k1sj;
          d1sj = dIsj; k1sj = i;
        } else if (dIsj>d2sj) {
          d2sj = dIsj; k2sj = i;
        } nIsj += 1.;
      }

      hJetNsj->Fill(dJet, nIsj, dNorm);
      if (d1sj>0.) { hJet1sj->Fill(dJet, d1sj, dNorm); hJet1sz->Fill(dJet, d1sj/dJet, dNorm); }
      if (d2sj>0.) { hJet2sj->Fill(dJet, d2sj, dNorm); hJet2sz->Fill(dJet, d2sj/dJet, dNorm); }

      if ((d1sj>0.) && (d2sj>0.)) {
        TVector3 v1sj; v1sj.SetPtEtaPhi(d1sj, trimmdSj[k1sj].eta(), trimmdSj[k1sj].phi());
        TVector3 v2sj; v2sj.SetPtEtaPhi(d2sj, trimmdSj[k2sj].eta(), trimmdSj[k2sj].phi());

        double dsj = d1sj - d2sj;
        double dsz = dsj / dJet;
        double dsr = v1sj.DeltaR(v2sj) / 2. / dJetR;

        hJetDsj->Fill(dJet, dsj, dNorm);
        hJetDsz->Fill(dJet, dsz, dNorm);
        hJetDsr->Fill(dJet, dsz, dsr, dNorm);
      }
    }
//=============================================================================

    delete evt;
    ascii_in >> evt;
  }
//=============================================================================

  TFile *file = TFile::Open(Form("%s.root",sFile.Data()), "NEW");
  list->Write();
  file->Close();
//=============================================================================

  cout << "DONE" << endl;
  return 0;
}

//.........这里部分代码省略.........
	  continue;
	
	TGraphErrors *g = etaGraphs[idxCent][idxDataSet][idxMoment];
	
	if (idxCent == 0) 
	  ShiftGraphX(g, -0.015);
	else if (idxCent == 1) 
	  ShiftGraphX(g, -0.005);
	else if (idxCent == 4) 
	  ShiftGraphX(g, 0.005);
	else if (idxCent == 8) 
	  ShiftGraphX(g, 0.015);
	
	ConfigGraph(g, idxMoment, idxCent);
	
	if (idxCent == 0) {
	  g->Draw("AP");
	  
	  if (idxMoment == 5) {
	    // TLine *line0 = new TLine(aMinX, 0, aMaxX, 0);
	    // line0->SetLineColor(kGray+1);
	    // line0->SetLineStyle(2);
	    // line0->SetLineWidth(2);
	    // line0->Draw();
	  }
	  else if (idxMoment == 6) {
	    TLine *line1 = new TLine(aMinX, 1, aMaxX, 1);
	    line1->SetLineColor(kGray+1);
	    line1->SetLineStyle(2);
	    line1->SetLineWidth(2);
	    line1->Draw();
	    legRat[idxDataSet]->AddEntry(line1, "Poisson", "l");
	  }
	}
	
	g->Draw("PSAME");
	
	if (idxMoment == 4) 
	  legRat[idxDataSet]->AddEntry(etaGraphs[idxCent][idxDataSet][4], Form("%s", cent1[idxCent]), "pl");
	
      } // for (int idxCent = 0; idxCent < 9; idxCent++) { 
    } // for (int idxMoment = 0 ; idxMoment < nMoments; ++idxMoment) {
    
    pad->cd(2);
    TLatex *texb_3 = new TLatex(0.05, 0.55, "Au+Au collisions #sqrt{#it{s}_{NN}} = 14.5 GeV");
    texb_3->SetTextSize(0.07);
    texb_3->Draw("same");
    
    TLatex *texb_3a = new TLatex(0.05, 0.49, "Net-Charge, 0.2 < #it{p}_{T} (GeV/#it{c}) < 2.0, 0-5%");
    texb_3a->SetTextSize(0.07);
    texb_3a->Draw("same");
    
    TLatex *texb_3b = new TLatex(0.05,0.44, "statistical errors only");
    texb_3b->SetTextSize(0.06);
    texb_3b->Draw("same");
    
    pad->cd(1);
    TLatex *texb_4 = new TLatex(0.7, 19, "STAR Preliminary");
    texb_4->SetTextSize(0.07);
    texb_4->Draw("same");
    
    pad->cd(3);
    legRat[idxDataSet]->Draw();
    
    pad->Modified();
  } // for (int idxDataSet = 0 ; idxDataSet < nDataSets; ++idxDataSet) {
  
  // -----------------------------------------------------

  SaveCanvas(name);

  // -----------------------------------------------------
  
  TFile *fOut = TFile::Open("STAR_QM2015_Preliminary.root", "UPDATE");
  fOut->cd();
  
  TList* list = new TList;

  for (int idxMoment = 4; idxMoment < nMoments; ++idxMoment) {
    for (int idxCent = 0; idxCent < nCent; ++idxCent) {
      if (idxCent > 1) 
	continue;

      if (idxCent == 0) 
	ShiftGraphX(etaGraphs[idxCent][0][idxMoment], 0.015);
      else if (idxCent == 1) 
	ShiftGraphX(etaGraphs[idxCent][0][idxMoment], 0.005);
      else if (idxCent == 4) 
	ShiftGraphX(etaGraphs[idxCent][0][idxMoment], -0.005);
      else if (idxCent == 8) 
	ShiftGraphX(etaGraphs[idxCent][0][idxMoment], -0.015);
      
      etaGraphs[idxCent][0][idxMoment]->SetName(Form("Net-Charge_%s_DeltaEta_14.5GeV_%s_stat", aMoments[idxMoment], cent[idxCent]));
      list->Add(etaGraphs[idxCent][0][idxMoment]);
    }
  }

  list->Write("Net-Charge_VsDeltaEta", TObject::kSingleKey);
  fOut->Close();
}

//.........这里部分代码省略.........
  std::vector<fastjet::PseudoJet> fjInput;
  const double dMass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
//=============================================================================

  TList *list = new TList();
  TH1D *hJet = new TH1D("hJet", "", 1000, 0., 1000.); hJet->Sumw2(); list->Add(hJet);

  enum                    {  kJet,   kMje,   kDsz,   kIsm,   kZsm,   kDsr, kVar };
  const TString sHist[] = { "aJet", "aMje", "aDsz", "aIsm", "aZsm", "aDsr" };
  const Int_t      nv[] = {   1000,    150,  120,      150,    150,    500 };
  const Double_t dMin[] = {     0.,     0.,   0.,       0.,     0.,     0. };
  const Double_t dMax[] = {  1000.,   150.,  1.2,     150.,    1.5,     5. };

  THnSparseD *hs = new THnSparseD("hs", "", kVar, nv, dMin, dMax); hs->Sumw2();
  for (Int_t i=0; i<kVar; i++) hs->GetAxis(i)->SetName(sHist[i].Data()); list->Add(hs);
//=============================================================================

  HepMC::IO_GenEvent ascii_in(Form("%s/%s.hepmc",sPath.Data(),sFile.Data()), std::ios::in);
  HepMC::GenEvent *evt = ascii_in.read_next_event();

  while (evt) {
    fjInput.resize(0);

    TLorentzVector vPar;
    for (HepMC::GenEvent::particle_const_iterator p=evt->particles_begin(); p!=evt->particles_end(); ++p) if ((*p)->status()==1) {
      vPar.SetPtEtaPhiM((*p)->momentum().perp(), (*p)->momentum().eta(), (*p)->momentum().phi(), dMass);

      if ((TMath::Abs(vPar.Eta())<dCutEtaMax)) {
        fjInput.push_back(fastjet::PseudoJet(vPar.Px(), vPar.Py(), vPar.Pz(), vPar.E()));
      }
    }
//=============================================================================

    fastjet::ClusterSequenceArea clustSeq(fjInput, jetsDef, areaDef);
    std::vector<fastjet::PseudoJet> includJets = clustSeq.inclusive_jets(dJetsPtMin);
//  std::vector<fastjet::PseudoJet> subtedJets = bkgSubtractor(includJets);
    std::vector<fastjet::PseudoJet> selectJets = selectJet(includJets);
//  std::vector<fastjet::PseudoJet> sortedJets = fastjet::sorted_by_pt(selectJets);

    for (int j=0; j<selectJets.size(); j++) {
      double dJet = selectJets[j].pt();

      hJet->Fill(dJet);
//=============================================================================

      fastjet::Filter trimmer(subjDef, fastjet::SelectorPtFractionMin(0.));
      fastjet::PseudoJet trimmdJet = trimmer(selectJets[j]);
      std::vector<fastjet::PseudoJet> trimmdSj = trimmdJet.pieces();

      double d1sj = -1.; int k1sj = -1;
      double d2sj = -1.; int k2sj = -1;
      for (int i=0; i<trimmdSj.size(); i++) {
        double dIsj = trimmdSj[i].pt(); if (dIsj<0.001) continue;

        if (dIsj>d1sj) {
          d2sj = d1sj; k2sj = k1sj;
          d1sj = dIsj; k1sj = i;
        } else if (dIsj>d2sj) {
          d2sj = dIsj; k2sj = i;
        }
      }

      if ((d1sj>0.) && (d2sj>0.)) {
        TLorentzVector v1sj; v1sj.SetPtEtaPhiM(d1sj, trimmdSj[k1sj].eta(), trimmdSj[k1sj].phi(), trimmdSj[k1sj].m());
        TLorentzVector v2sj; v2sj.SetPtEtaPhiM(d2sj, trimmdSj[k2sj].eta(), trimmdSj[k2sj].phi(), trimmdSj[k2sj].m());
        TLorentzVector vIsj = v1sj + v2sj;

        Double_t dIsm = vIsj.M();
        Double_t dMje = selectJets[j].m();
        Double_t dVar[] = { dJet, dMje, (d1sj-d2sj)/dJet, dIsm, dIsm/dMje, v1sj.DeltaR(v2sj)/2./dJetR };

        hs->Fill(dVar);
      }
    }
//=============================================================================

    delete evt;
    ascii_in >> evt;
  }
//=============================================================================

  TString sXsec = sFile; sXsec.ReplaceAll("out", "xsecs");
  TFile *file = TFile::Open(Form("%s/xsecs/%s.root",sPath.Data(),sXsec.Data()), "READ");
  TH1D *hPtHat        = (TH1D*)file->Get("hPtHat");        hPtHat->SetDirectory(0);
  TH1D *hWeightSum    = (TH1D*)file->Get("hWeightSum");    hWeightSum->SetDirectory(0);
  TProfile *hSigmaGen = (TProfile*)file->Get("hSigmaGen"); hSigmaGen->SetDirectory(0);
  file->Close();
//=============================================================================

  file = TFile::Open(Form("%s.root",sFile.Data()), "NEW");
  hPtHat->Write();
  hWeightSum->Write();
  hSigmaGen->Write();
  list->Write();
  file->Close();
//=============================================================================

  cout << "DONE" << endl;
  return 0;
}

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

TCanvas* cZPC_Mean_Uncalib = new TCanvas("cZPC_Mean_Uncalib","cZPC_Mean_Uncalib", 0,0,750,900);
hzpcUncalib->Draw("ep");
cZPC_Mean_Uncalib->Print(Form("%s/cZPC_Mean_Uncalib.png",plotDir.Data()));

TCanvas* cZPA_Mean_Uncalib = new TCanvas("cZPA_Mean_Uncalib","cZPA_Mean_Uncalib", 0,0,750,900);
hzpaUncalib->Draw("ep");
cZPA_Mean_Uncalib->Print(Form("%s/cZPA_Mean_Uncalib.png",plotDir.Data()));

TCanvas* cZEM1_Mean_Values = new TCanvas("cZEM1_Mean_Values","cZEM1_Mean_Values", 0,0,750,900);
hzem1->Draw("ep");
cZEM1_Mean_Values->Print(Form("%s/cZEM1_Mean_Values.png",plotDir.Data()));

TCanvas* cZEM2_Mean_Values = new TCanvas("cZEM2_Mean_Values","cZEM2_Mean_Values", 0,0,750,900);
hzem2->Draw("ep");
cZEM2_Mean_Values->Print(Form("%s/cZEM2_Mean_Values.png",plotDir.Data()));

//---------------------------------------------------------------------------------------------------
//centroids
//---------------------------------------------------------------------------------------------------
TCanvas* cZNA_X_centroid = new TCanvas("cZNA_X_centroid","cZNA_X_centroid", 0,0,750,900);
hzna_Xcentroid->Draw();
cZNA_X_centroid->Print(Form("%s/cZNA_X_centroid.png",plotDir.Data()));

TCanvas* cZNA_Y_centroid = new TCanvas("cZNA_Y_centroid","cZNA_Y_centroid", 0,0,750,900);
hzna_Ycentroid->Draw();
cZNA_Y_centroid->Print(Form("%s/cZNA_Y_centroid.png",plotDir.Data()));

TCanvas* cZNC_X_centroid = new TCanvas("cZNC_X_centroid","cZNC_X_centroid", 0,0,750,900);
hznc_Xcentroid->Draw();
cZNC_X_centroid->Print(Form("%s/cZNC_X_centroid.png",plotDir.Data()));

TCanvas* cZNC_Y_centroid = new TCanvas("cZNC_Y_centroid","cZNC_Y_centroid", 0,0,750,900);
hznc_Ycentroid->Draw();
cZNC_Y_centroid->Print(Form("%s/cZNC_Y_centroid.png",plotDir.Data()));

//---------------------------------------------------------------------------------
//timing
//---------------------------------------------------------------------------------
TCanvas* cTimingSum = new TCanvas("cTimingSum","cTimingSum",0,0,750,900);
hzn_TDC_Sum->Draw();
cTimingSum->Print(Form("%s/cTimingSum.png",plotDir.Data()));

TCanvas* cTimingDiff = new TCanvas("cTimingDiff","cTimingDiff",0,0,750,900);
hzn_TDC_Diff->Draw();
cTimingDiff->Print(Form("%s/cTimingDiff.png",plotDir.Data()));

//----------------------------------------------------------------------
//out
//----------------------------------------------------------------------
printf(" preparing output tree\n");
tree->Branch("run",&runNumber,"runNumber/I");
tree->Branch("ZNC_mean_value",&ZNC_mean,"ZNC_mean/D");
tree->Branch("ZNA_mean_value",&ZNA_mean,"ZNA_mean/D");
tree->Branch("ZPC_mean_value",&ZPC_mean,"ZPC_mean/D");
tree->Branch("ZPA_mean_value",&ZPA_mean,"ZPA_mean/D");
tree->Branch("ZNC_mean_uncalib",&ZNCuncalib_mean,"ZNCuncalib_mean/D");
tree->Branch("ZNA_mean_uncalib",&ZNAuncalib_mean,"ZNAuncalib_mean/D");
tree->Branch("ZPC_mean_uncalib",&ZPCuncalib_mean,"ZPCuncalib_mean/D");
tree->Branch("ZPA_mean_uncalib",&ZPAuncalib_mean,"ZPAuncalib_mean/D");
tree->Branch("ZEM1_mean_value",&ZEM1_mean,"ZEM1_mean/D");
tree->Branch("ZEM2_mean_value",&ZEM2_mean,"ZEM2_mean/D");
tree->Branch("ZNC_X_Centroid",&ZNC_XCent,"ZNC_XCent/D");
tree->Branch("ZNC_Y_Centroid",&ZNC_YCent,"ZNC_YCent/D");
tree->Branch("ZNA_X_Centroid",&ZNA_XCent,"ZNA_XCent/D");
tree->Branch("ZNA_Y_Centroid",&ZNA_YCent,"ZNA_YCent/D");
tree->Branch("ZNC_X_Centroid_Err",&ZNC_XCent_err,"ZNC_XCent_err/D");
tree->Branch("ZNC_Y_Centroid_Err",&ZNC_YCent_err,"ZNC_YCent_err/D");
tree->Branch("ZNA_X_Centroid_Err",&ZNA_XCent_err,"ZNA_XCent_err/D");
tree->Branch("ZNA_Y_Centroid_Err",&ZNA_YCent_err,"ZNA_YCent_err/D");
tree->Branch("ZN_TDC_Sum",&ZN_TDC_Sum,"ZN_TDC_Sum/D");
tree->Branch("ZN_TDC_Diff",&ZN_TDC_Diff,"ZN_TDC_Diff/D");
tree->Branch("ZN_TDC_Sum_Err",&ZN_TDC_Sum_err,"ZN_TDC_Sum_err/D");
tree->Branch("ZN_TDC_Diff_Err",&ZN_TDC_Diff_err,"ZN_TDC_Diff_err/D");
tree->Fill();

if(hZNCpmcUncalib) list.Add(cZNC_Spectra_Uncal);
if(hZNApmcUncalib) list.Add(cZNA_Spectra_Uncal);
if(hZPCpmcUncalib) list.Add(cZPC_Spectra_Uncal);
if(hZPApmcUncalib) list.Add(cZPA_Spectra_Uncal);
list.Add(cZEM1_Spectra);
list.Add(cZEM2_Spectra);
list.Add(cTimingSum);
list.Add(cTimingDiff);
list.Add(cZNA_X_centroid);
list.Add(cZNA_Y_centroid);
list.Add(cZNC_X_centroid);
list.Add(cZNC_Y_centroid);

TFile *fout;
fout = TFile::Open("prodQAhistos.root", "update");
if(!fout) fout = new TFile("prodQAhistos.root");
fout->cd();
list.Write();
tree->Write();
fout->Close();

return 0;

}

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

    hT0T0CVsRun->SetBinContent(irun+1,StartTime_pT0CT0);
    hT0T0CVsRun->SetBinError(irun+1,StartTime_pT0CT0Err);
    hT0T0CVsRun->GetXaxis()->SetBinLabel(irun+1,runlabel);
    
    hT0BestVsRunRes->SetBinContent(irun+1,StartTime_pBestT0_Res);
    hT0BestVsRunRes->SetBinError(irun+1, 1.e-5);
    hT0BestVsRunRes->GetXaxis()->SetBinLabel(irun+1,runlabel);

    hT0FillVsRunRes->SetBinContent(irun+1,StartTime_pFillT0_Res);
    hT0FillVsRunRes->SetBinError(irun+1, 1.e-5);   
    hT0FillVsRunRes->GetXaxis()->SetBinLabel(irun+1,runlabel);

    hT0TOFVsRunRes->SetBinContent(irun+1,StartTime_pTOFT0_Res);
    hT0TOFVsRunRes->SetBinError(irun+1, 1.e-5);
    hT0TOFVsRunRes->GetXaxis()->SetBinLabel(irun+1,runlabel);

    hT0T0ACVsRunRes->SetBinContent(irun+1,StartTime_pT0ACT0_Res);
    hT0T0ACVsRunRes->SetBinError(irun+1, 1.e-5);
    hT0T0ACVsRunRes->GetXaxis()->SetBinLabel(irun+1,runlabel);

    hT0T0AVsRunRes->SetBinContent(irun+1,StartTime_pT0AT0_Res);
    hT0T0AVsRunRes->SetBinError(irun+1, 1.e-5);
    hT0T0AVsRunRes->GetXaxis()->SetBinLabel(irun+1,runlabel);

    hT0T0CVsRunRes->SetBinContent(irun+1,StartTime_pT0CT0_Res);
    hT0T0CVsRunRes->SetBinError(irun+1, 1.e-5);
    hT0T0CVsRunRes->GetXaxis()->SetBinLabel(irun+1,runlabel);

  }
  
  TFile * fout=new TFile(outfilename,"recreate");
  fout->cd();
  lista.Write();
  fout->Close();
    
  gStyle->SetOptStat(10);

  TString plotext = "png";
  const TString desiredext = gSystem->Getenv("TOFQAPLOTEXTENSION");
  if(desiredext.EqualTo("pdf") || desiredext.EqualTo("root")) plotext = desiredext; 
  else if(!desiredext.IsNull()) cout<<"Unrecognized extension: '"<<desiredext<<"'"<<endl;
  
  //Plot t-texp trend
  TCanvas* cPeakDiffTimeVsRun = new TCanvas("cPeakDiffTimeVsRun","cPeakDiffTimeVsRun", 50,50,1050, 550);
  hPeakDiffTimeVsRun->GetYaxis()->SetRangeUser(-50.,50.);
  hPeakDiffTimeVsRun->Draw();
  cPeakDiffTimeVsRun->Print(Form("%s/cPeakDiffTimeVsRun.%s", plotDir.Data(), plotext.Data()));
	
  TCanvas* cSpreadDiffTimeVsRun = new TCanvas("cSpreadDiffTimeVsRun","cSpreadDiffTimeVsRun", 50,50,1050, 550);
  hSpreadDiffTimeVsRun->GetYaxis()->SetRangeUser(0.,400.);
  hSpreadDiffTimeVsRun->Draw();
  cSpreadDiffTimeVsRun->Print(Form("%s/cSpreadDiffTimeVsRun.%s", plotDir.Data(), plotext.Data()));

  //Plot average of t-texp-t0tof and resolution trend
  TCanvas* cMeanTOFResVsRun = new TCanvas("cMeanTOFResVsRun","cMeanTOFResVsRun", 50,50,1050, 550);
  hMeanTOFResVsRun->GetYaxis()->SetRangeUser(-50.,50.);
  hMeanTOFResVsRun->Draw();
  cMeanTOFResVsRun->Print(Form("%s/cMeanTOFResVsRun.%s", plotDir.Data(), plotext.Data()));
       
  TCanvas* cSigmaTOFResVsRun = new TCanvas("cSigmaTOFResVsRun","cSigmaTOFResVsRun", 50,50,1050, 550);
  hSigmaTOFResVsRun->GetYaxis()->SetRangeUser(0.,200.);
  hSigmaTOFResVsRun->Draw();
  cSigmaTOFResVsRun->Print(Form("%s/cSigmaTOFResVsRun.%s", plotDir.Data(), plotext.Data()));

  //Plot matching efficiency trend

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

    AliGenPythiaEventHeader *pHeadPy = (AliGenPythiaEventHeader*)pHeader->GenEventHeader();

    if (!pHeadPy) continue;
    hPtHat->Fill(pHeadPy->GetPtHard());
//=============================================================================

    for (Int_t i=0; i<pStack->GetNtrack(); i++) if (pStack->IsPhysicalPrimary(i)) {
      TParticle *pTrk = pStack->Particle(i); if (!pTrk) continue;
      if (TMath::Abs(pTrk->Eta())>dCutEtaMax) { pTrk = 0; continue; }
//    TParticlePDG *pPDG = pTrk->GetPDG(); if (!pPDG) { pTrk = 0; continue; }

      fjInput.push_back(fastjet::PseudoJet(pTrk->Px(), pTrk->Py(), pTrk->Pz(), pTrk->P()));

//    pPDG = 0;
      pTrk = 0;
    }
//=============================================================================

    fastjet::ClusterSequenceArea clustSeq(fjInput, jetsDef, areaDef);
    std::vector<fastjet::PseudoJet> includJets = clustSeq.inclusive_jets(dJetsPtMin);
//  std::vector<fastjet::PseudoJet> subtedJets = bkgSubtractor(includJets);
    std::vector<fastjet::PseudoJet> selectJets = selectJet(includJets);
//  std::vector<fastjet::PseudoJet> sortedJets = fastjet::sorted_by_pt(selectJets);

    for (int j=0; j<selectJets.size(); j++) {
      double dJet = selectJets[j].pt();

      hJet->Fill(dJet);
//=============================================================================

      fastjet::Filter trimmer(subjDef, fastjet::SelectorPtFractionMin(0.));
      fastjet::PseudoJet trimmdJet = trimmer(selectJets[j]);
      std::vector<fastjet::PseudoJet> trimmdSj = trimmdJet.pieces();

      double nIsj = 0.;
      double d1sj = -1.; int k1sj = -1;
      double d2sj = -1.; int k2sj = -1;
      for (int i=0; i<trimmdSj.size(); i++) {
        double dIsj = trimmdSj[i].pt(); if (dIsj<0.001) continue;

        hJetIsj->Fill(dJet, dIsj);
        hJetIsz->Fill(dJet, dIsj/dJet);

        if (dIsj>d1sj) {
          d2sj = d1sj; k2sj = k1sj;
          d1sj = dIsj; k1sj = i;
        } else if (dIsj>d2sj) {
          d2sj = dIsj; k2sj = i;
        } nIsj += 1.;
      }

      hJetNsj->Fill(dJet, nIsj);
      if (d1sj>0.) { hJet1sj->Fill(dJet, d1sj); hJet1sz->Fill(dJet, d1sj/dJet); }
      if (d2sj>0.) { hJet2sj->Fill(dJet, d2sj); hJet2sz->Fill(dJet, d2sj/dJet); }

      if ((d1sj>0.) && (d2sj>0.)) {
        double dsj = d1sj - d2sj;
        double dsz = dsj / dJet;

        hJetDsj->Fill(dJet, dsj);
        hJetDsz->Fill(dJet, dsz);
      }
    }
//=============================================================================

    pStack  = 0;
    pHeadPy = 0;
    pHeader = 0;
  }
//=============================================================================

  rl->UnloadgAlice();
  rl->UnloadHeader();
  rl->UnloadKinematics();
  rl->RemoveEventFolder();
//=============================================================================

  TFile *file = TFile::Open(Form("%s/pyxsec_hists.root",sPath.Data()), "READ");
  TList *lXsc = (TList*)file->Get("cFilterList");
  file->Close();

  TH1D     *hWeightSum = (TH1D*)lXsc->FindObject("h1Trials");   hWeightSum->SetName("hWeightSum");
  TProfile *hSigmaGen  = (TProfile*)lXsc->FindObject("h1Xsec"); hSigmaGen->SetName("hSigmaGen");
//=============================================================================

  file = TFile::Open(Form("%s.root",sFile.Data()), "NEW");
  hPtHat->Write();
  hWeightSum->Write();
  hSigmaGen->Write();
  list->Write();
  file->Close();
//=============================================================================

  cout << "DONE" << endl;
//=============================================================================

  return 0;
}

//.........这里部分代码省略.........
      fHistZvtxRes2[2]->Draw("sames");
      leg2->Draw();
      ccr->SaveAs("plotXYZVtxRMS.png");

      TCanvas *ccc = new TCanvas("ccc", "ccc", 1200, 800);
      ccc->Divide(3, 2);
      ccc->cd(1);
      fHistDecLen[0]->GetYaxis()->SetTitleOffset(1.45);
      fHistDecLen[0]->Draw();
      fHistDecLen[2]->Draw("sames");
      if (drawOnlyDzerDplus == 0)
	fHistDecLen[4]->Draw("sames");
      if (drawOnlyDzerDplus == 0)
	fHistDecLen[6]->Draw("sames");
      if (drawOnlyDzerDplus == 0)
	fHistDecLen[8]->Draw("sames");
      leg3->Draw();

      ccc->cd(2);
      fHistCosp[0]->GetYaxis()->SetTitleOffset(1.45);
      fHistCosp[0]->Draw();
      fHistCosp[2]->Draw("sames");
      if (drawOnlyDzerDplus == 0)
	fHistCosp[4]->Draw("sames");
      if (drawOnlyDzerDplus == 0)
	fHistCosp[6]->Draw("sames");
      if (drawOnlyDzerDplus == 0)
	fHistCosp[8]->Draw("sames");
      leg4->Draw();

      ccc->cd(3);
      hEffPt[0]->GetYaxis()->SetTitleOffset(1.45);
      hEffPt[0]->Draw();
      hEffPt[2]->Draw("sames");
      if (drawOnlyDzerDplus == 0)
	hEffPt[4]->Draw("sames");
      if (drawOnlyDzerDplus == 0)
	hEffPt[6]->Draw("sames");
      if (drawOnlyDzerDplus == 0)
	hEffPt[8]->Draw("sames");
      leg3->Draw();

      ccc->cd(4);
      fHistDecLen[1]->GetYaxis()->SetTitleOffset(1.45);
      fHistDecLen[1]->Draw();
      fHistDecLen[3]->Draw("sames");
      if (drawOnlyDzerDplus == 0)
	fHistDecLen[5]->Draw("sames");
      if (drawOnlyDzerDplus == 0)
	fHistDecLen[7]->Draw("sames");
      if (drawOnlyDzerDplus == 0)
	fHistDecLen[9]->Draw("sames");
      leg3->Draw();

      ccc->cd(5);
      fHistCosp[1]->GetYaxis()->SetTitleOffset(1.45);
      fHistCosp[1]->Draw();
      fHistCosp[3]->Draw("sames");
      if (drawOnlyDzerDplus == 0)
	fHistCosp[5]->Draw("sames");
      if (drawOnlyDzerDplus == 0)
	fHistCosp[7]->Draw("sames");
      if (drawOnlyDzerDplus == 0)
	fHistCosp[9]->Draw("sames");
      leg4->Draw();

      ccc->cd(6);
      hEffPt[1]->GetYaxis()->SetTitleOffset(1.45);
      hEffPt[1]->Draw();
      hEffPt[3]->Draw("sames");
      if (drawOnlyDzerDplus == 0)
	hEffPt[5]->Draw("sames");
      if (drawOnlyDzerDplus == 0)
	hEffPt[7]->Draw("sames");
      if (drawOnlyDzerDplus == 0)
	hEffPt[9]->Draw("sames");
      leg3->Draw();

      ccc->SaveAs("plot_DL_cosp_Eff_prompt_fd.png");
    }
  }

  trtree->Fill();

  if(runNumber>0){
    TFile* foutfile=new TFile("trendingHF.root","recreate");
    trtree->Write();
    TDirectory* outdir=foutfile->mkdir(dirD2H->GetName());
    outdir->cd();
    listD2H->Write(listD2H->GetName(),1);
    foutfile->cd();
    if(dir && list){
      TDirectory* outdir2=foutfile->mkdir(dir->GetName());
      outdir2->cd();
      list->Write(list->GetName(),1);
    }
    foutfile->Close();
    delete foutfile;
  }
}