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

void treegraph(TString filename) {
   gROOT->SetStyle("Plain");
   gStyle->SetOptDate();


   Double_t x, y;
   Int_t nlines = 0;
   TFile *f = new TFile("graph.root","RECREATE");

   TCanvas *canvas_graph = new TCanvas("canvas_graph", "y vs x",467,89,400,700);

   TTree t;
   t.ReadFile(filename,"x:y");
   t.Draw("x:y","","goff");

   TGraph *g = new TGraph(t.GetSelectedRows(),t.GetV1(),t.GetV2());
   g->SetTitle(filename+": Y vs X");
   g->GetXaxis()->SetTitle("x[a.u.]");
   g->GetYaxis()->SetTitle("y[a.u.]");

   g->SetMarkerStyle(21);
   g->Draw("AP");

   f->Write();
}

void CheckOutput(){
  //
  //
  //
  TFile * f  = TFile::Open("Filtered.root");
  TTree * highPt = (TTree*)f->Get("highPt");
  TTree * treeV0s = (TTree*)f->Get("V0s");
  //
  // Export variable:
  //
  Double_t ratioHighPtV0Entries=treeV0s->GetEntries()/Double_t(treeV0s->GetEntries()+highPt->GetEntries()+0.000001);
  Double_t ratioHighPtV0Size=treeV0s->GetZipBytes()/Double_t(treeV0s->GetZipBytes()+highPt->GetZipBytes()+0.000001);
  printf("#UnitTest:\tAliAnalysisTaskFiltered\tRatioPtToV0Entries\t%f\n",ratioHighPtV0Entries);
  printf("#UnitTest:\tAliAnalysisTaskFiltered\tRatioPtToV0Size\t%f\n",ratioHighPtV0Size);
  //
  //
  Double_t ratioPointsV0 = 2*treeV0s->GetBranch("friendTrack0.fCalibContainer")->GetZipBytes()/Double_t(0.00001+treeV0s->GetZipBytes());
  Double_t ratioPointsHighPt = highPt->GetBranch("friendTrack.fCalibContainer")->GetZipBytes()/Double_t(0.00001+highPt->GetZipBytes());
  printf("#UnitTest:\tAliAnalysisTaskFiltered\tRatioPointsV0\t%f\n",ratioPointsV0);
  printf("#UnitTest:\tAliAnalysisTaskFiltered\tRatioPointsHighPt\t%f\n",ratioPointsHighPt);
  //
  // a.) Check track correspondence
  //
  Int_t entries= highPt->Draw("(friendTrack.fTPCOut.fP[3]-esdTrack.fIp.fP[3])/sqrt(friendTrack.fTPCOut.fC[9]+esdTrack.fIp.fC[9])","friendTrack.fTPCOut.fP[3]!=0","");
  // here we should check if the tracks
  Double_t pulls=TMath::RMS(entries, highPt->GetV1());
  printf("#UnitTest:\tAliAnalysisTaskFiltered\tFriendPull\t%2.4f\n",pulls);
  printf("#UnitTest:\tAliAnalysisTaskFiltered\tFriendOK\t%d\n",pulls<10);

}

/*
 * The Ferromagnet_Scan function reads in
 * a data file of magentic field scans
 * from within the ferromagnet and returns
 * a graph of the magnetic permeability
 * of the ferromagnet vs the external field
 * provided by the Helmholtz coil.
 */
TGraphErrors* plot_uvB(
			 const TString scan_file,
			 TF1* calib_fit,
			 double R,
			 double R_sig
)
{

  /*Read in Data File to ROOT Tree*/
  cout<< "processing file " << scan_file <<endl;
  TTree *TData = new TTree();
  TData->ReadFile(scan_file, "t/D:I:B");
  /*Use TTree Draw command to write branches to usable arrays*/
  int n = TData->Draw("I:TMath::Abs(B):0.0:0.005", "", "goff");

  vector<double> B_ext, B_in, Bratio, u;
  for(int i = 0; i < n; i++)
    {
      B_ext.push_back( calib_fit->Eval(TData->GetV1()[i]) );
      B_in.push_back( TData->GetV2()[i] );
      Bratio.push_back( B_in[i] / B_ext[i]);
      u.push_back( (Bratio[i]*(R**2) + Bratio[i] - 2 - 2*sqrt((Bratio[i]**2)*(R**2) - Bratio[i]*(R**2) - Bratio[i] + 1 ) ) / (Bratio[i]*(R**2) - Bratio[i]) );

    }

  TGraphErrors *g_uvB = new TGraphErrors(n, &B_ext[0], &u[0], TData->GetV3(), TData->GetV4());

  /*  
  g_uvB->Fit("pol1", "", "", 10, 60);
  cout << "Permeability at 50: " << g_uvB->GetFunction("pol1")->Eval(50) << endl;
  cout << "Permeability at 40: " << g_uvB->GetFunction("pol1")->Eval(40) << endl;
  cout << "Permeability at 30: " << g_uvB->GetFunction("pol1")->Eval(30) << endl;
  */
  return g_uvB;
}

/*!
 * ### Example usage:
\code
    .L $AliRoot_SRC/STAT/Macros/aliExternalInfo.C+
    drawLogbook("LHC15o","pass1","QA.TPC.meanTPCncl>0","runDuration:totalEventsPhysics:totalNumberOfFilesMigrated:MonALISA.RCT.tpc_value");
    drawLogbook("LHC10h","pass2","totalEventsPhysics>1000&&totalNumberOfFilesMigrated>20","runDuration:totalEventsPhysics:totalNumberOfFilesMigrated:MonALISA.RCT.tpc_value");
    drawLogbook("LHC11h","pass2","totalEventsPhysics>1000&&totalNumberOfFilesMigrated>20","runDuration:totalEventsPhysics:totalNumberOfFilesMigrated:MonALISA.RCT.tpc_value");
\endcode
*/
void drawLogbook(TString period, TString pass,TString runSelection="", TString varSelection="runDuration:totalEventsPhysics:totalNumberOfFilesMigrated"){
  TTree * treeLogbook = info.GetTree("Logbook",period,pass,"QA.TPC;QA.TRD;QA.ITS;MonALISA.RCT");
  TObjArray *varList=0;
  if (varSelection[0]=='['){ //variable list using class selection
    // Use class selection to select variables
    varList=AliTreePlayer::selectMetadata(treeLogbook,varSelection,0);
    Int_t nvars=varList->GetEntries();
    for (Int_t i=0; i<nvars;i++){
      TString vname=varList->At(i)->GetName();
      vname.ReplaceAll(".class","");
    }
    // varList=AliTreePlayer::selectMetadata(treeLogbook, "[class==\"Logbook&&Time\"]",0);
  }else{
    varList=varSelection.Tokenize(":");
  }
  Int_t nvars=varList->GetEntries();
  TCanvas * canvas = new TCanvas("drawLogbook","drawLogbook",1200,1000);
  canvas->Divide(1,nvars);
  for (Int_t i=0; i<nvars; i++){
    canvas->cd(i+1);
    TString var=TString::Format("%s:run",varList->At(i)->GetName());
    TStatToolkit::MakeGraphSparse(treeLogbook,var,runSelection,25,1,1)->Draw("ap");
    Double_t mean = TMath::Mean(treeLogbook->GetSelectedRows(),treeLogbook->GetV1());
    Double_t sum  = treeLogbook->GetSelectedRows()*mean;
    latex.DrawLatexNDC(0.15,0.8,TString::Format("Mean=%0.0f",mean));
    latex.DrawLatexNDC(0.15,0.7,TString::Format("Sum=%0.0f",sum));
  }
  canvas->SaveAs(TString::Format("%s_%s.png",pass.Data(),period.Data()).Data());
}

void histogram() {
	// Read data from ASCII file and create histogram/ntuple combo
	ifstream in;
 	in.open("/users/ronnie/git/Journal-Analysis/data.txt");

	Text_t month;
 	Float_t date, day, words;
	TFile *file = new TFile("histogramTest.root", "CREATE");
	TH1F *histo = new TH1F("histo", "writing distribution", 100, 0, 2);
	TTree *Tree = new TTree("ntuple", "data from file");
					
	Tree->ReadFile(Form("/users/ronnie/git/Journal-Analysis/Daily Journals/data.txt"), "month:date:day:words");
	graph = new TGraph(200, &words, &day);
	// Get the total number of words written and print out to CLI
	Tree->Print();
	Int_t numEntries = (Int_t)(Tree->GetEntries());
	TBranch *wordsBranch = Tree->GetBranch("words");
	for(Int_t i = 0; i < numEntries; i++) {
		Tree->GetEntry(i, 0);
	}
	// Draw the Tree (words vs. day) on canvas
	Tree->SetEstimate(Tree->GetEntries());
	Tree->Draw("words:day");
	// Have to loop back through drawn tree to get values and count number of words total
	Float_t numWords = 0.0;
	Double_t *array = Tree->GetV1();
	for (Int_t i = 0; i < numEntries; i++) {
		numWords += array[i];
	}
	cout << numWords << endl;
}

makePlot_check_all()
{

  gStyle->SetOptStat(0);

  vector< TString > v_infile;
  //  v_infile.push_back( TString("output/fitResults_run001_1layer_sheath.root") );
  //  v_infile.push_back( TString("output/fitResults_run002_1layer_sheath.root") );
  v_infile.push_back( TString("output/fitResults_run003_1layer_sheath.root") );
  v_infile.push_back( TString("output/fitResults_run004_2layer_sheath.root") );
  v_infile.push_back( TString("output/fitResults_run005_2layer_sheath.root") );
  v_infile.push_back( TString("output/fitResults_run006_3layer_sheath.root") );
  v_infile.push_back( TString("output/fitResults_run007_4layer_sheath.root") );
  v_infile.push_back( TString("output/fitResults_run008_2layer_cylinder.root") );
  v_infile.push_back( TString("output/fitResults_run009_2layer_cylinder.root") );
  //  v_infile.push_back( TString("output/fitResults_run010_1layer_sheath.root") );
  v_infile.push_back( TString("output/fitResults_run011_1layer_sheath.root") );
  v_infile.push_back( TString("output/fitResults_run012_1layer_sheath.root") );
  v_infile.push_back( TString("output/fitResults_run013_1layer_sheath.root") );
  v_infile.push_back( TString("output/fitResults_run014_1layer_sheath.root") );
  v_infile.push_back( TString("output/fitResults_run015_1layer_sheath.root") );
  v_infile.push_back( TString("output/fitResults_run016_1layer_sheath.root") );
  v_infile.push_back( TString("output/fitResults_run017_2layer_sheath.root") );
  v_infile.push_back( TString("output/fitResults_run018_1layer_sheath.root") );
  v_infile.push_back( TString("output/fitResults_run019_1layer_sheath.root") );
  v_infile.push_back( TString("output/fitResults_run020_1layer_sheath.root") );
  v_infile.push_back( TString("output/fitResults_run021_1layer.root") );
  for ( unsigned f = 0; f < v_infile.size(); f++ )
    {
      cout << "Plot " << v_infile.at(f) << endl;

      TCanvas *c0 = new TCanvas();

      TFile *fin = new TFile( v_infile.at(f) );
      TTree *tin = (TTree*)fin->Get("fitResults");
      tin->Draw("bshield:bext:bshield_err:bext_err");

      TGraphErrors *gshield = new TGraphErrors( tin->GetEntries(), tin->GetV2(), tin->GetV1(), tin->GetV4(), tin->GetV3() );
      gshield->SetTitle(v_infile.at(f));
      gshield->GetXaxis()->SetTitle("B_{ext} [mT]");
      gshield->GetYaxis()->SetTitle("B_{shield} [mT]");

      gshield->Draw("AP");

      fin->Close();

      if ( f == 0 )
        c0->Print("plot_check_all.ps(");
      else if ( f == v_infile.size() - 1 )
        c0->Print("plot_check_all.ps)");
      else
        c0->Print("plot_check_all.ps");

    }



}

/// \brief Cache  MC production trees, store summary information in formated text files -> root trees
/// \param dataType  -
/// \param fileList
void CacheTestMCProductions(TString dataType, const char *fileList=NULL){
  AliExternalInfo info;
  info.fLoadMetadata=kFALSE;
  TObjArray* periodList = NULL;
  TArrayI nRuns;
  if (fileList!=NULL) {
    periodList=(gSystem->GetFromPipe(TString::Format("cat %s", fileList).Data())).Tokenize("\n");
    nRuns.Set(periodList->GetEntries());

  }else{
    TTree * tree = info.GetTree("MonALISA.ProductionMC","","");
    Int_t nProd=tree->GetEntries();
    periodList = new TObjArray(nProd);
    nRuns.Set(nProd);
    TLeaf *leaf = tree->GetLeaf("Tag");
    TLeaf *leafRuns = tree->GetLeaf("Number_of_runs");
    for (Int_t iProd=0; iProd<nProd; iProd++){
      tree->GetEntry(iProd);
      TString prodName=((char*)leaf->GetValuePointer());
      if (prodName.Contains("LHC")==0) continue;
      periodList->AddAt(new TObjString(((char*)leaf->GetValuePointer())),iProd);
      nRuns[iProd]=leafRuns->GetValue();
    }
    delete tree;
  }
  for (Int_t iPeriod=0; iPeriod<periodList->GetEntriesFast(); iPeriod++){
    TObjString * pName= (TObjString*)periodList->At(iPeriod);
    if (pName==NULL) continue;
    TTree* tree = info.GetTree(dataType.Data(),periodList->At(iPeriod)->GetName(),"passMC");
    if (tree){
      Int_t entries=tree->Draw("run","1","goff");
      TString sInfo=periodList->At(iPeriod)->GetName();
      sInfo+="\t";
      sInfo+=dataType;
      sInfo+="\t";
      sInfo+=TString::Format("%d\t",entries);
      sInfo+=TString::Format("%d\t",nRuns[iPeriod]);
      for (Int_t j=0; j<entries; j++) {
        sInfo+=TString::Format("%2.0f,",tree->GetV1()[j]);
        ::Info("CacheTestMCProductionsRun:","%s\t%s\t%d\t%d\t%d\t%2.0f",periodList->At(iPeriod)->GetName(),dataType.Data(),entries,nRuns[iPeriod],j, tree->GetV1()[j]);
      }
      sInfo+="0";
      ::Info("CacheTestMCProductionsPeriod:","%s\n",sInfo.Data());
      delete tree;
    }else{
      ::Error("CacheTestMCProductionsPeriod:","%s\t%s\t-1\t%d\t0",periodList->At(iPeriod)->GetName(), dataType.Data(),nRuns[iPeriod]);
    }
  }
}

/*
 * The Calibration function reads in the
 * calibration file and returns a graph
 * that is used to find the relation
 * between current and magentic field
 * for the Helmholtz coil.
 */
TF1* Calibrate(
		    const char* f_calib
		    )
{
  /*Read in Calibration File*/
  cout<< "processing file " << f_calib <<endl;
  TTree *TCalib = new TTree();
  TCalib->ReadFile(f_calib, "t/D:I:B");
  int n = TCalib->Draw("I:TMath::Abs(B)", "", "goff");
  TGraph *g_calib = new TGraph(n, TCalib->GetV1(), TCalib->GetV2() );
  g_calib->SetTitle("");
  // g_calib->Draw("AP");
  g_calib->Fit("pol1", "q");
  TF1* calib_fit = g_calib->GetFunction("pol1");

  return calib_fit;

}

TGraphErrors* plot_ramp(
			 const TString scan_file,
			 double R,
			 double R_sig
)
{

  /*Read in Data File to ROOT Tree*/
  cout<< "processing file " << scan_file <<endl;
  TTree *TData = new TTree();
  TData->ReadFile(scan_file, "t/D:x:y:z:B1:B1Range:B2:B2Range:B3:B3Range:V:T1:T2:T3:T4:T5:T6:T7:T8");
  /*Use TTree Draw command to write branches to usable vector*/
  int n = TData->Draw("TMath::Abs(B2):TMath::Abs(B3):0.0:0.005", "", "goff");

  vector<double> B_ext, B_in, Bratio, u, u_err;
  for(int i = 0; i < n; i++)
    {
  //     if(i == 0)
  //     {
  //     	B_ext.push_back(100.00);
  //     }
  //     else{
  //     	  B_ext.push_back( TData->GetV1()[i]);
  //     }
      B_ext.push_back( TData->GetV1()[i] );
      B_in.push_back( TData->GetV2()[i]);
      // B_in.push_back( 90.53);
      // if(TMath::Abs(B_in[i]) <= TMath::Abs(B_ext[i]))
      // {
	Bratio.push_back( B_in[i] / B_ext[i]);
	u.push_back( (Bratio[i]*(R**2) + Bratio[i] - 2 - 2*sqrt((Bratio[i]**2)*(R**2) - Bratio[i]*(R**2) - Bratio[i] + 1 ) ) / (Bratio[i]*(R**2) - Bratio[i]) );
	// u_err.push_back( u[i]*( 5.5*((0.005/Bratio[i])**2.) + 10*((R_sig/R)**2.)  )**0.5 );
	u_err.push_back( 0.01 );
	// u_err.push_back( u[i]*( 5.5*((0.01/Bratio[i])**2.))  )**0.5 );
	// }
      }

  TGraphErrors *g_uvB = new TGraphErrors(n, &B_ext[0], &u[0], TData->GetV3(), &u_err[0]);

  return g_uvB;
}

int evaluate( std::string filelist, std::string outfile )
{
  gStyle->SetOptStat(0);

  TCanvas *ctemp = new TCanvas();

  TCanvas *cres = new TCanvas("TimeDependence");
  TH1F* hres = new TH1F("hres","",100,0,650);
  hres->GetYaxis()->SetRangeUser(0,50);
  hres->SetTitle("");
  hres->GetXaxis()->SetTitle("time (s)");
  hres->GetYaxis()->SetTitle("B_{int} (mT)");
  hres->Draw();
  leg = new TLegend(0.2,0.6,0.9,0.9);
//  leg->SetHeader("The Legend Title"); // option "C" allows to center the header
  leg->SetNColumns(5);

  vector< double > v_Bint;
  vector< double > v_BintErr;
  vector< double > v_Bext;
  vector< double > v_BextErr;

  /* Loop over all lines in input file */
  std::ifstream infilelist(filelist);
  std::string line;

  unsigned colorcounter=38;

  while (std::getline(infilelist, line))
  {
    // skip lines with '#' and empty lines
    if ( line.find("#") != string::npos )
      {
        cout << "Skip line " << line << endl;
        continue;
      }

    if ( line == "" )
      continue;

    //cout << "Processing file " << line << endl;


    TString infilename("data_calib/");
    infilename.Append(line);

    TFile *fin = new TFile( infilename );
    TTree *tin = (TTree*)fin->Get("t");

    ctemp->cd();
    tin->Draw("Bi:time");
    TGraph *gtime = new TGraph(tin->GetEntries(), &(tin->GetV2()[0]), &(tin->GetV1()[0]));
    gtime->SetLineColor(colorcounter);
    colorcounter++;

    TH1F* hBext = new TH1F("hBext","",100,0,1000);
    tin->Draw("Bo >> hBext");

    cres->cd();
    gtime->Draw("lsame");

    double Bext_i = hBext->GetMean();
    double BextErr_i = hBext->GetRMS();

    double Bint_i = gtime->Eval(590);
    double BintErr_i = 0;

    /* add legend entry */
    TString legname("B_ext ~ ");
    legname += (int)Bext_i;
    leg->AddEntry(gtime,legname,"l");

    cout << "B_ext: " << Bext_i << " \t B_int: " << Bint_i << endl;

    v_Bint.push_back(Bint_i);
    v_BintErr.push_back(BintErr_i);
    v_Bext.push_back(Bext_i);
    v_BextErr.push_back(BextErr_i);

  }

  cres->cd();
  leg->Draw();

  TGraphErrors *gfinal = new TGraphErrors(v_Bext.size(), &(v_Bext[0]), &(v_Bint[0]), &(v_BextErr[0]), &(v_BintErr[0]));
  gfinal->Sort();
  gfinal->SetName("Bint_Vs_Bext");
  gfinal->SetTitle("");
  gfinal->GetXaxis()->SetTitle("B_{ext} (mT)");
  gfinal->GetYaxis()->SetTitle("B_{int} (mT)");

  TCanvas *cfinal = new TCanvas();
  gfinal->Draw("APL");

  /* Save output graph */
  TString outfilename("output/");
  outfilename.Append(outfile);
  TFile *fout = new TFile(outfilename,"RECREATE");

  cres->Write();
//.........这里部分代码省略.........

void fillGraphFromTreeVar(std::string& treedrawspec,int index,wGraph_t *&pwg)
{
  // Sample treedrawspec:
  // mytree:"TDCwinstart[%d]:runnum","evtnum==1","P"
  //
  vector<string> v_tokens;
  string tid;
  TString drawspec;
  Tokenize(treedrawspec,v_tokens,":",true);
  if( (v_tokens.size() < 2) ||
      (!v_tokens[0].size())  ||
      (!v_tokens[2].size())    ) {
    cerr << "malformed root tree draw spec treeid:\"varexp\",\"selection\",option: " << treedrawspec << endl;
    return;
  }

  tid = v_tokens[0];
  for (size_t i=2; i<v_tokens.size(); i++) {
    drawspec += v_tokens[i];
  }
  int fmtcnt = drawspec.CountChar('%');

  if (fmtcnt) { // use index for tree array var
    switch(fmtcnt) {
    case 1: drawspec = Form(drawspec,index); break;
    case 2: drawspec = Form(drawspec,index,index); break;
    case 3: drawspec = Form(drawspec,index,index,index); break;
    case 4: drawspec = Form(drawspec,index,index,index,index); break;
    case 5: drawspec = Form(drawspec,index,index,index,index,index); break;
    case 6: drawspec = Form(drawspec,index,index,index,index,index,index); break;
    default:
      cerr << "More than six fmt specifiers in drawspec found, fix me! " << drawspec <<endl;
      exit(-1);
    }
  }
  if( gl_verbose)
    cout<<"drawspec="<<drawspec<<endl;

  TTree *tree = findTree(tid);
  assert (tree);

  // can't use comma as delimiter since histo with binning spec may be supplied
  TObjArray *tokens = drawspec.Tokenize("\"");
  TString hname;

  TString varexp = ((TObjString *)(*tokens)[0])->GetString();

  if( gl_verbose)
    cout<<"varexp="<<varexp<<endl;
  switch(tokens->GetEntriesFast()) {
  case 1:
    tree->Draw(varexp,"","goff");
    break;
  case 3:
    {
      TString cut = ((TObjString *)(*tokens)[2])->GetString();
      tree->Draw(varexp,cut,"goff"); 
    }
    break;
  case 4: // assume the cut string is blank
    {
      TString gopt = ((TObjString *)(*tokens)[3])->GetString();
      gopt = gopt + " goff";
      tree->Draw(varexp,"",gopt);
    }
    break;
  case 5:
    {
      TString cut  = ((TObjString *)(*tokens)[2])->GetString();
      TString gopt = ((TObjString *)(*tokens)[4])->GetString();
      gopt = gopt + " goff";
      tree->Draw(varexp,cut,gopt);
    }
    break;
  default:
    cerr << "malformed root tree draw spec treeid:varexp,selection,option";
    for (int i=0; i<tokens->GetEntriesFast(); i++)
      cerr << i<<": "<< ((TObjString *)(*tokens)[i])->GetString() << " ";
    cerr << endl;
    break;
  }

  assert(tree->GetSelectedRows());

  if (!pwg)
    pwg = new wGraph_t();

  assert(pwg);
  pwg->gr = new TGraph(tree->GetSelectedRows(),
		       tree->GetV2(), tree->GetV1());

}                                                // fillGraphFromTreeVar

/*#include "TTree.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TMultiGraph.h"
#include "TRint.h"
#include <time.h>
#include <stdio.h>

const short MaxNN = 6;
*/
void make_scatterplot(TString data_file, TString plot_name, TString header)
{
	//TApplication program = new TRInt(); 
	//pull in data
	TTree *t = new TTree(); 
	t->ReadFile(data_file);
	t->SetName("t");
	
	//TString name = "run2.png";
	TCanvas *BG = new TCanvas("c1", "Read Velocity on Local Disk for CMS3 Files", 1920, 1080);
	BG->cd();
	TPad *c = new TPad("MainPad", "My main pad", 0, 0, 1, 0.9);
	c->Draw();
	c->Divide(2,2);

	//In first slot have Time vs. Buffer for 1 Concurrent Read
	c->cd(1);
	int n1 = t->Draw("VelocityMBps:BufferSize", "ConcurrentReads==1", "goff");
	TGraph *ghist1 = new TGraph(n1, t->GetV2(), t->GetV1());
	ghist1->SetName("ghist1");
	ghist1->SetMarkerStyle(3);
	ghist1->SetMarkerColor(1);
	ghist1->SetTitle("Single File Read");
	ghist1->GetXaxis()->SetTitle("Buffer Size (Bytes)");
	ghist1->GetYaxis()->SetTitle("Read Velocity (MB/s)");
	ghist1->Draw("ap");
	
	//In second slot have Time vs. Buffer for 3 Concurrent Reads
	c->cd(2);
	int n2 = t->Draw("VelocityMBps:BufferSize", "ConcurrentReads==3", "goff");
	TGraph *ghist2 = new TGraph(n2, t->GetV2(), t->GetV1());
	Double_t DP2x[n2], DP2y[n2];
	ghist2->SetName("ghist2");
	ghist2->SetMarkerStyle(3);
	ghist2->SetMarkerColor(1);
	ghist2->SetTitle("3 Concurrent Reads");
	ghist2->GetXaxis()->SetTitle("Buffer Size (Bytes)");
	ghist2->GetYaxis()->SetTitle("Read Velocity (MB/s)");
	ghist2->Draw("ap");
	
	//In third slot have Time vs. Buffer for 3 Concurrent Reads
	c->cd(3);
	int n3 = t->Draw("VelocityMBps:BufferSize", "ConcurrentReads==6", "goff");
	TGraph *ghist3 = new TGraph(n3, t->GetV2(), t->GetV1());
	Double_t DP3x[n3], DP3y[n3];
	ghist3->SetName("ghist3");
	ghist3->SetMarkerStyle(3);
	ghist3->SetMarkerColor(1);
	ghist3->SetTitle("6 Concurrent Reads");
	ghist3->GetXaxis()->SetTitle("Buffer Size (Bytes)");
	ghist3->GetYaxis()->SetTitle("Read Velocity (MB/s)");
	ghist3->Draw("ap");
	
	//In fourth slot have Time vs. Buffer for 10 Concurrent Reads
	c->cd(4);
	int n4 = t->Draw("VelocityMBps:BufferSize", "ConcurrentReads==10", "goff");
	TGraph *ghist4 = new TGraph(n4, t->GetV2(), t->GetV1());
	ghist4->SetMarkerStyle(3);
	ghist4->SetName("ghist4");
	ghist4->SetMarkerColor(1);
	ghist4->SetTitle("10 Concurrent Reads");
	ghist4->GetXaxis()->SetTitle("Buffer Size (Bytes)");
	ghist4->GetYaxis()->SetTitle("Read Velocity (MB/s)");
	ghist4->Draw("ap");
	
	//Draw to screen
	//Initialize Canvas
	
	c->cd(0);
	c->Draw();
	
	BG->cd();
	TText *title = new TText(.5,.95, header);
	title->SetTextAlign(22);
	title->Draw();
	
	gDirectory->Add(ghist1);
	gDirectory->Add(ghist2);
	gDirectory->Add(ghist3);
	gDirectory->Add(ghist4);
	gDirectory->Add(t);
	gPad->SaveAs(plot_name);
}

/// Cache MC production information
void CacheTrendingProductions(TString dataType){
  AliExternalInfo info;
  info.fLoadMetadata=kFALSE;
  TObjArray* periodList = NULL, *idList=NULL;
  //
  TTree * tree = info.GetTree("MonALISA.ProductionCycle","","");
  Int_t nProd=tree->GetEntries();
  periodList = new TObjArray(nProd);
  idList= new TObjArray(nProd);
  TLeaf *leafTag = tree->GetLeaf("Tag");
  TLeaf *leafID   =  tree->GetLeaf("ID");
  for (Int_t iProd=0; iProd<nProd; iProd++){
    tree->GetEntry(iProd);
    TString prodName=((char*)leafTag->GetValuePointer());
    TString  idName =TString::Format("%d",TMath::Nint(leafID->GetValue()));
    if (prodName.Contains("LHC")==0) continue;
    periodList->AddAt(new TObjString(prodName),iProd);
    idList->AddAt(new TObjString(idName),iProd);
  }
  delete tree;
  //
  for (Int_t iPeriod=0; iPeriod<periodList->GetEntriesFast(); iPeriod++) {
    TObjString * pName= (TObjString*)idList->At(iPeriod);
    if (pName==NULL) continue;
    TTree* treeP = info.GetTreeProdCycleByID(idList->At(iPeriod)->GetName());
    if (treeP==NULL) continue;
    TLeaf *leafOutput = treeP->GetLeaf("outputdir");
    Int_t nRuns= treeP->GetEntries();
    treeP->GetEntry(0);
    TString path=((char*)leafOutput->GetValuePointer());
    TObjArray *pArray = path.Tokenize("/");
    if (pArray==NULL) continue;
    Int_t nElems=pArray->GetEntries();
    if (nElems<4) continue;
    TString aperiod=pArray->At(3)->GetName();
    TString apass =pArray->At(nElems-1)->GetName();
    delete pArray;
    ::Info("CacheTrendingProductions","%s\t%s\t%s\t%s\t%d",idList->At(iPeriod)->GetName(),path.Data(), aperiod.Data(),apass.Data(),nRuns);
    delete treeP;
    TTree* treeQA = info.GetTree(dataType.Data(),aperiod.Data(),apass.Data());
    if (treeQA){
      Int_t entries=treeQA->Draw("run","1","goff");
      TString sInfo=aperiod;
      sInfo+="\t";
      sInfo+=apass;
      sInfo+="\t";
      sInfo+=dataType;
      sInfo+="\t";
      sInfo+=TString::Format("%d\t",entries);
      sInfo+=TString::Format("%d\t",nRuns);
      for (Int_t j=0; j<entries; j++) {
        sInfo+=TString::Format("%2.0f,",treeQA->GetV1()[j]);
        ::Info("CacheTrendingProductionsRun:","%s\t%s\t%s\t%d\t%d\t%2.0f",aperiod.Data(),apass.Data(),dataType.Data(),entries,nRuns,treeQA->GetV1()[j]);
      }
      sInfo+="0";
      ::Info("CacheTrendingProductionsPeriod:","%s\n",sInfo.Data());
      delete treeQA;
    }else{
      ::Error("CacheTrendingProductionsPeriod:","%s\t%s\t%s\t-1\t%d\t0",aperiod.Data(),apass.Data(), dataType.Data(),nRuns);
    }
  }
}

int
plot_compare_fun4all_eicroot()
{
  /**** Chose input files ****/
  /* File with IR mangets configuration*/
  TString fname_irmag("example/proton-magnets-250GeV-opt2.dat");
  /* File with Fun4All output*/
  TString fname_fun4all("example/eRHIC_proton-magnets-250GeV-opt2_250GeV_0mrad.root");
  /* File with EICROOT output*/
  TString fname_eicroot("example/eicroot-track_proton-magnets-250GeV-opt2_250GeV_0mrad.txt");

  /* Open iput file with trajectories from GEANT4 */
  TFile *fin = new TFile(fname_fun4all);

  /* Get tree from file */
  TTree *tin = (TTree*)fin->Get("T");

  int nhits = 0;
  tin->SetBranchAddress("n_G4HIT_FWDDISC",&nhits);

  /* create graph of particle trajectory */
  /* Use only first event (for now) */
  tin->GetEntry(0);
  cout << "hits: " << nhits << endl;

  tin->Draw("G4HIT_FWDDISC.x:G4HIT_FWDDISC.z","Entry$==0","");

  TGraph* g1 = new TGraph(nhits*2, &(tin->GetV2()[0]), &(tin->GetV1()[0]));
  g1->SetMarkerStyle(7);
  g1->SetMarkerSize(1);
  g1->SetMarkerColor(kRed);

  /* Get tree from file */
  TTree *tin2 = new TTree();
  tin2->ReadFile(fname_eicroot,"x/F:y:z");

  int nhits = tin2->GetEntries();
  tin2->Draw("x:z","","");

  TGraph* g2 = new TGraph(nhits, &(tin2->GetV2()[0]), &(tin2->GetV1()[0]));
  g2->SetMarkerStyle(7);
  g2->SetMarkerSize(1);
  g2->SetMarkerColor(kGreen+1);

  /* Create frame histogram for plot */
  TH1F *h1 = new TH1F("h1","",10,0,15000);
  h1->GetXaxis()->SetRangeUser(0,5000);
  h1->GetYaxis()->SetRangeUser(-50,70);
  h1->GetXaxis()->SetTitle("Z(cm)");
  h1->GetYaxis()->SetTitle("X(cm)");

  /* Plot frame histogram */
  TCanvas *c1 = new TCanvas();
  h1->Draw("AXIS");

  /* Read IR configuration file- this needs to go somewhere else using parameters and a .root file to store them */
  ifstream irstream(fname_irmag);

  if(!irstream.is_open())
    {
      cout << "ERROR: Could not open IR configuration file " << fname_irmag << endl;
      return -1;
    }

  while(!irstream.eof()){
    string str;
    getline(irstream, str);
    if(str[0] == '#') continue; //for comments
    stringstream ss(str);

    string name;
    double center_z, center_x, center_y, aperture_radius, length, angle, B, gradient;

    ss >> name >> center_z >> center_x >> center_y >> aperture_radius >> length >> angle >> B >> gradient;

    if ( name == "" ) continue; //for empty lines

    // convert units from m to cm
    center_x *= 100;
    center_y *= 100;
    center_z *= 100;
    aperture_radius *= 100;
    length *= 100;

    // define magnet outer radius
    float outer_radius = 30.0; // cm

    //flip sign of dipole field component- positive y axis in Geant4 is defined as 'up',
    // positive z axis  as the hadron-going direction
    // in a right-handed coordinate system x,y,z
    B *= -1;

    // convert angle from millirad to rad
    angle = (angle / 1000.);

    // Place IR component
    cout << "New IR component: " << name << " at z = " << center_z << endl;

    string volname = "IRMAGNET_";
    volname.append(name);
//.........这里部分代码省略.........

int
TrackParametrization( TString csvfile="fitslices_out.csv" )
{

  /* Read data from input file */
  TTree *tres = new TTree();
  tres->ReadFile( csvfile, "ptrue:etatrue:psig:psig_err:pmean:pmean_err:norm", ',' );

  /* Print read-in tree */
  tres->Print();

  /* colors array */
  unsigned colors[8] = {1,2,3,4,6,7,14,16};

  /* Create vector of theta values to include for visualization*/
  vector< double > etas_vis;
  etas_vis.push_back(-2.75);
  etas_vis.push_back(-2.25);
  etas_vis.push_back(-1.75);
  etas_vis.push_back(-0.25);
  etas_vis.push_back( 0.25);
  etas_vis.push_back( 1.75);
  etas_vis.push_back( 2.25);

//  etas_vis.push_back(-3.25);
//  etas_vis.push_back(-2.25);
//  etas_vis.push_back(-1.25);
//  etas_vis.push_back(-0.25);
//  etas_vis.push_back( 0.25);
//  etas_vis.push_back( 1.25);
//  etas_vis.push_back( 2.25);
//  etas_vis.push_back( 3.25);

  /* Create vector of theta values to include for fitting*/
  vector< double > etas_fit;
  for ( double eta = -4.45; eta < 4.5; eta += 0.1 )
    etas_fit.push_back( eta );

  /* Create fit function */
  TF1* f_momres = new TF1("f_momres", "sqrt( [0]*[0] + [1]*[1]*x*x )" );

  cout << "\nFit function: " << f_momres->GetTitle() << "\n" << endl;

  /* Create scratch canvas */
  TCanvas *cscratch = new TCanvas("cscratch");

  /* Create framehistogram */
  TH1F* hframe = new TH1F("hframe","",100,0,40);
  hframe->GetYaxis()->SetRangeUser(0,0.15);
  hframe->GetYaxis()->SetNdivisions(505);
  hframe->GetXaxis()->SetTitle("Momentum (GeV/c)");
  hframe->GetYaxis()->SetTitle("#sigma_{p}/p");

  /* create combined canvas plot */
  TCanvas *c1 = new TCanvas();
  hframe->Draw();

  /* Create legend */
  TLegend* leg_eta = new TLegend( 0.2, 0.6, 0.5, 0.9);
  leg_eta->SetNColumns(2);

  /* Create ofstream to write fit parameter results */
  ofstream ofsfit("track_momres_new.csv");
  ofsfit<<"eta,par1,par1err,par2,par2err"<<endl;

  /* Create resolution-vs-momentum plot with fits for each selected theta value */
  for ( int i = 0; i < etas_fit.size(); i++ )
    {
      /* Switch to scratch canvas */
      cscratch->cd();

      double eta = etas_fit.at(i);

      /* No tracking outside -4 < eta < 4 */
      if ( eta < -4 || eta > 4 )
	continue;

      cout << "\n***Eta = " << eta << endl;

      /* Define range of theta because float comparison with fixed value doesn't work
	 too well for cuts in ROOT trees */
      double eta_min = eta * 0.999;
      double eta_max = eta * 1.001;

      /* Cut for tree */
      TCut cutx( Form("ptrue > 1 && ( (etatrue > 0 && (etatrue > %f && etatrue < %f)) || (etatrue < 0 && (etatrue < %f && etatrue > %f)) )", eta_min, eta_max, eta_min, eta_max) );

      /* "Draw" tree on scratch canvas to fill V1...V4 arrays */
      tres->Draw("psig:ptrue:psig_err:0", cutx );

      /* Create TGraphErrors with selected data from tree */
      TGraphErrors *gres = new TGraphErrors( tres->GetEntries(cutx),
					     &(tres->GetV2())[0],
					     &(tres->GetV1())[0],
					     &(tres->GetV4())[0],
					     &(tres->GetV3())[0] );

      /* reset function parameters before fit */
      f_momres->SetParameter(0,0.1);
      f_momres->SetParameter(1,0.1);
//.........这里部分代码省略.........