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

void PlotParsDistr(TFile* f, TTree* tr, TString strMillepedeRes, TString strOutdir)
{

  for (int isubd=PXB; isubd<=TEC; isubd++){

    TString canvName="c_";
    canvName+=strMillepedeRes;
    canvName+="_";
    canvName+=StrPlotType(PARS);
    canvName+="_";
    canvName+=subdLabels[isubd];
    canvName.ReplaceAll(".res","");

    f->cd();
    TCanvas* canv = new TCanvas(canvName,canvName,600,600);
    canv->Divide(3,3);

    for (int parInd=1; parInd<=9; parInd++){
      canv->cd(parInd);
      TString strCut="((label%20-1)%9+1)==";
      strCut+=parInd;
      strCut+=" && label<700000 && ";
      strCut+=StrCutSubd(isubd);
      TString hName="hPars_";
      hName+=subdLabels[isubd];
      hName+="_";
      hName+= StrPar(parInd);
      
      TTree* trCut = tr->CopyTree(strCut);
      float up = trCut->GetMaximum("parVal");
      float low = trCut->GetMinimum("parVal");
      std::cout<<"low="<<low<<", up="<<up<<", nent="<<trCut->GetEntries()<<std::endl;
      TH1F* h = new TH1F(hName,hName,100,10000*low,10000*up);
      TString strDraw="10000*parVal>>";
      strDraw+=hName;
      trCut->Draw(strDraw,strCut,"goff");
      h->SetMarkerStyle(2);
      h->Draw("EP");
    }// end of loop over parInd

    canvName+=".png";
    TString saveName=strOutdir+canvName;
    canv->SaveAs(saveName);
    saveName.ReplaceAll(".png",".pdf");
    canv->SaveAs(saveName);

  }//end of loop over isubd

}// end of PlotParsDistr

TGraphErrors* PlotLightYieldGraph()
{
  string filename;
  vector<double> x,y,ex,ey;
  while(1){
    cout<<"\n\nEnter next file to process; <enter> to finish."<<endl;
    getline(cin, filename);
    if(filename=="")
      break;
    
    //load the tree
    TTree* Events = GetEventsTree(filename.c_str());
    if(!Events)
      continue;
    gROOT->ProcessLine(".x analysis/Aliases.C");
    double start = Events->GetMinimum("timestamp");
    double end = Events->GetMaximum("timestamp");
    double error = 0;
    TString title = TString("Na22 Spectrum, ") + 
      TString(filename)(TRegexp("Run......"));
    TH1F* hist = new TH1F("Na22Spec",title,200,1000,2500);
    Events->Draw("sumspec >> Na22Spec","min > 0","e");
    double yield = Fit511Photopeak(hist,&error);
    
    x.push_back((start+end)/2.);
    ex.push_back((end-start)/2.);
    y.push_back(yield);
    ey.push_back(error);
  }
  if(x.size() == 0){
    cerr<<"No valid points found!"<<endl;
    return 0;
  }
  TGraphErrors* graph = new TGraphErrors(x.size(),&x[0],&y[0],&ex[0],&ey[0]);
  graph->Draw("ape");
  TAxis* xax = graph->GetXaxis();
  xax->SetTitle("Run Time");
  xax->SetTimeDisplay(1);
  xax->SetTimeFormat("%Y/%m/%d");
  xax->SetTimeOffset(1,"gmt");
  TAxis* yax = graph->GetYaxis();
  yax->SetTitle("511 keV Light Yield [pe/keV]");
  graph->Draw("ape");
  return graph;
}

void nbx_check_2(const char * filename="counts.root")
{
  TFile * infile = new TFile(filename,"READ");
  TTree * tr = (TTree*) infile->Get("sca");

  Int_t NRUNS_tmp = tr->GetMaximum("i");
  const Int_t NRUNS = NRUNS_tmp;
  TH1D * h[NRUNS];
  char h_n[NRUNS][256];
  char h_t[NRUNS][256];
  char cut[NRUNS][256];

  TCanvas * c = new TCanvas("c","c",1400,1000);

  Double_t max;

  for(Int_t i=0; i<NRUNS; i++)
  {
    printf("i=%d\n",i);

    sprintf(h_n[i],"nbx_bx_%d",i+1);
    sprintf(h_t[i],"N_{bx} vs. bXing for i=%d",i+1);
    h[i] = new TH1D(h_n[i],h_n[i],120,0,120);
    sprintf(cut[i],"tot_bx*(i==%d)",i+1);
    tr->Project(h_n[i],"bx",cut[i]);

    max = h[i]->GetMaximum();
    h[i]->Scale(1/max);

    c->Clear();
    c->SetGrid(1,1);
    h[i]->Draw();

    if(i==0) c->Print("nbx_check/nbx_vs_bxing.pdf(","pdf");
    else if(i+1==NRUNS) c->Print("nbx_check/nbx_vs_bxing.pdf)","pdf");
    else c->Print("nbx_check/nbx_vs_bxing.pdf");
  };
};

void paracoor( TString fin = "TMVA.root", Bool_t useTMVAStyle = kTRUE )
{
   // set style and remove existing canvas'
   TMVAGlob::Initialize( useTMVAStyle );

   // checks if file with name "fin" is already open, and if not opens one
   TFile* file = TMVAGlob::OpenFile( fin );  
   TTree* tree = (TTree*)file->Get("TestTree");
   if(!tree) {
      cout << "--- No TestTree saved in ROOT file. Parallel coordinates will not be plotted" << endl;
      return;
   }

   // first get list of leaves in tree
   TObjArray* leafList = tree->GetListOfLeaves();
   vector<TString> vars;
   vector<TString> mvas;   
   for (Int_t iar=0; iar<leafList->GetSize(); iar++) {
      TLeaf* leaf = (TLeaf*)leafList->At(iar);
      if (leaf != 0) {
         TString leafName = leaf->GetName();
         if (leafName != "type" && leafName != "weight"  && leafName != "boostweight" &&
             leafName != "class" && leafName != "className" && leafName != "classID" && 
             !leafName.Contains("prob_")) {
            // is MVA ?
            if (TMVAGlob::ExistMethodName( leafName )) {
               mvas.push_back( leafName );
            }
            else {
               vars.push_back( leafName );
            }
         }
      }
   }

   cout << "--- Found: " << vars.size() << " variables" << endl;
   cout << "--- Found: " << mvas.size() << " MVA(s)" << endl;
   

   TString type[2] = { "Signal", "Background" };
   const Int_t nmva = mvas.size();
   TCanvas* csig[nmva];
   TCanvas* cbkg[nmva];
   for (Int_t imva=0; imva<mvas.size(); imva++) {
      cout << "--- Plotting parallel coordinates for : " << mvas[imva] << " & input variables" << endl;

      for (Int_t itype=0; itype<2; itype++) {

         // create draw option
         TString varstr = mvas[imva] + ":";
         for (Int_t ivar=0; ivar<vars.size(); ivar++) varstr += vars[ivar] + ":";
         varstr.Resize( varstr.Last( ':' ) );

         // create canvas
         TString mvashort = mvas[imva]; mvashort.ReplaceAll("MVA_","");
         TCanvas* c1 = (itype == 0) ? csig[imva] : cbkg[imva];
         c1 = new TCanvas( Form( "c1_%i",itype ), 
                           Form( "Parallel coordinate representation for %s and input variables (%s events)", 
                                 mvashort.Data(), type[itype].Data() ), 
                           50*(itype), 50*(itype), 750, 500 );      
         tree->Draw( varstr.Data(), Form("classID==%i",1-itype) , "para" );
         c1->ToggleEditor();
         gStyle->SetOptTitle(0);

         TParallelCoord*    para   = (TParallelCoord*)gPad->GetListOfPrimitives()->FindObject( "ParaCoord" );
         TParallelCoordVar* mvavar = (TParallelCoordVar*)para->GetVarList()->FindObject( mvas[imva] );
         Double_t minrange = tree->GetMinimum( mvavar->GetName() );
         Double_t maxrange = tree->GetMaximum( mvavar->GetName() );
         Double_t width    = 0.2*(maxrange - minrange);
         Double_t x1 = minrange, x2 = x1 + width;
         TParallelCoordRange* parrange = new TParallelCoordRange( mvavar, x1, x2 );
         parrange->SetLineColor(4);
         mvavar->AddRange( parrange );

         para->AddSelection("-1");

         for (Int_t ivar=1; ivar<TMath::Min(Int_t(vars.size()) + 1,3); ivar++) {
            TParallelCoordVar* var = (TParallelCoordVar*)para->GetVarList()->FindObject( vars[ivar] );
            minrange = tree->GetMinimum( var->GetName() );
            maxrange = tree->GetMaximum( var->GetName() );
            width    = 0.2*(maxrange - minrange);

            switch (ivar) {
            case 0: { x1 = minrange; x2 = x1 + width; break; }
            case 1: { x1 = 0.5*(maxrange + minrange - width)*0.02; x2 = x1 + width*0.02; break; }
            case 2: { x1 = maxrange - width; x2 = x1 + width; break; }
            }

            parrange = new TParallelCoordRange( var, x1, x2 );
            parrange->SetLineColor( ivar == 0 ? 2 : ivar == 1 ? 5 : 6 );
            var->AddRange( parrange );

            para->AddSelection( Form("%i",ivar) );
         }

         c1->Update();

         TString fname = Form( "plots/paracoor_c%i_%s", imva, itype == 0 ? "S" : "B" );
         TMVAGlob::imgconv( c1, fname );
      }
//.........这里部分代码省略.........

//--------------------------------------
void FBI3D(TString namefile="", TString namefile2=""){

  Int_t    MaxNCoinc = 9e7;
  Int_t*   POS = (Int_t*)malloc(MaxNCoinc*sizeof(Int_t));
  FILE* fichero;
  FILE* fichero2;

  cout << namefile<<" "<<namefile2<<endl;

// === STEP 1 === NORMALIZATION ================

if (namefile2=="normf.raw"){   // Allows reusing previous normalization image
  cout << "Using precalculated Normalization file" << endl;
  fichero2 = fopen(namefile2,"rb");
  fread(SENS,sizeof(float),nvoxels,fichero2);
  fclose(fichero2); 
}else{          // Create normalization from ROOT file (filtered to reduce noise)  
  int nt = 10;  // Number of filter iterations
  SENS = normalization(namefile2,nt);
  // Storing the normalization image so that it can be used for other reconstructions of this scanner
  fichero2 = fopen("normf.raw","wb");    
  fwrite(SENS,sizeof(float),nvoxels,fichero2);
  fclose(fichero2);
}

// === STEP 2 === READ ROOT FILE  ======

 int posic; 
 TFile *f1 = new TFile(namefile,"READ");
 TTree *Coincidences = (TTree*)f1->Get("Coincidences"); 
 ncoinc = Coincidences->GetEntries();
 cout << "Simulated Coincidences = " << ncoinc << endl;
 define_branches(Coincidences);
 Scanner = (Coincidences->GetMaximum("globalPosX1")) - (Coincidences->GetMinimum("globalPosX1")); // Diameter of the scanner from the detected counts
 FOV_Z = (Coincidences->GetMaximum("globalPosZ1")) - (Coincidences->GetMinimum("globalPosZ1"));
 dz = float(NZS)/FOV_Z;  

// === STEP 3 === REFERENCE IMAGE (OPTIONAL) ====

for(Int_t iV=0; iV<nvoxels ; iV++){IMG_N[iV]=0.;}  // Initial Image
for(Int_t ic = 0; ic <= ncoinc; ic++){     // 
  Coincidences->GetEntry(ic);
  ix = int(source_X*dxy + RESM);
  iy = int(source_Y*dxy + RESM);
  iz = int(source_Z*dz + NZM);      
  if (ix>=0 && ix<RES && iy>=0 && iy<RES && iz>=0 && iz<NZS){   // Check voxel valid
   iV = iz*RES*RES+iy*RES+ix;
   IMG_N[iV]++;
  }
 }

fichero = fopen("image_ref.raw","wb");
fwrite(IMG_N,sizeof(float),nvoxels,fichero);
fclose(fichero);

Float_t sens_corr = 0.;
for(Int_t iV=0; iV<nvoxels ; iV++){
 sens_corr = SENS[iV];
 if (sens_corr>0){sens_corr = 1.0/sens_corr;}
 IMG_N[iV]*=sens_corr;
}  // Reference image corrected by sensitivity
fichero = fopen("image_ref_senscor.raw","wb");
fwrite(IMG_N,sizeof(float),nvoxels,fichero);
fclose(fichero);

// === STEP 4 === INITIAL IMAGE ================
for(Int_t iV=0; iV<nvoxels ; iV++){IMG[iV]=0;}  // Initial Image
for(Int_t iV=0; iV<ncoinc; iV++){POS[iV]=-1;}  // Initial Position

// === STEP 5 === ITERATIONS ==================
for(Int_t iter=1; iter<=niter ; iter++){        // Iterations (1 iteration consists of a whole loop over all data)
for(Int_t ic = 0;  ic<ncoinc;  ic++){     
 Coincidences->GetEntry(ic); 
 if ( ic % 100000 == 0 ) cout << "Iter= " << iter << " Processing " << ic << " coincidences " << endl;  
  posic = POS[ic]; 
  POS[ic] = Arek(posic,det1_X,det1_Y,det1_Z,det2_X,det2_Y,det2_Z,IMG,SENS);
 } // COINCIDENCES
} // ITERATIONS

f1->Close();

fichero = fopen("image.raw","wb");
if(fichero){
 fwrite(IMG,sizeof(int),nvoxels,fichero);
 fclose(fichero);
}

// === STEP 6 === FINAL SENSITIVY NORMALIZATION ====
for(Int_t iV=0; iV<nvoxels ; iV++){
  if (SENS[iV]>0){
    IMG_N[iV] = float(IMG[iV])/SENS[iV];
  }else{
    IMG_N[iV] = 0.;
  }
}  

fichero = fopen("image_norm.raw","wb");
if(fichero){
 fwrite(IMG_N,sizeof(float),nvoxels,fichero);
//.........这里部分代码省略.........

void PlotDecisionBoundary( TString weightFile = "weights/TMVAClassification_BDT.weights.xml",TString v0="var0", TString v1="var1", TString dataFileName = "/home/hvoss/TMVA/TMVA_data/data/data_circ.root") 
{   
   //---------------------------------------------------------------
   // default MVA methods to be trained + tested

   // this loads the library
   TMVA::Tools::Instance();

   std::cout << std::endl;
   std::cout << "==> Start TMVAClassificationApplication" << std::endl;


   //
   // create the Reader object
   //
   TMVA::Reader *reader = new TMVA::Reader( "!Color:!Silent" );    

   // create a set of variables and declare them to the reader
   // - the variable names must corresponds in name and type to 
   // those given in the weight file(s) that you use
   Float_t var0, var1;
   reader->AddVariable( v0,                &var0 );
   reader->AddVariable( v1,                &var1 );

   //
   // book the MVA method
   //
   reader->BookMVA( "M1", weightFile ); 
   
   TFile *f = new TFile(dataFileName);
   TTree *signal     = (TTree*)f->Get("TreeS");
   TTree *background = (TTree*)f->Get("TreeB");


//Declaration of leaves types
   Float_t         svar0;
   Float_t         svar1;
   Float_t         bvar0;
   Float_t         bvar1;
   Float_t         sWeight=1.0; // just in case you have weight defined, also set these branchaddresses
   Float_t         bWeight=1.0*signal->GetEntries()/background->GetEntries(); // just in case you have weight defined, also set these branchaddresses

   // Set branch addresses.
   signal->SetBranchAddress(v0,&svar0);
   signal->SetBranchAddress(v1,&svar1);
   background->SetBranchAddress(v0,&bvar0);
   background->SetBranchAddress(v1,&bvar1);


   UInt_t nbin = 50;
   Float_t xmax = signal->GetMaximum(v0.Data());
   Float_t xmin = signal->GetMinimum(v0.Data());
   Float_t ymax = signal->GetMaximum(v1.Data());
   Float_t ymin = signal->GetMinimum(v1.Data());
 
   xmax = TMath::Max(xmax,background->GetMaximum(v0.Data()));  
   xmin = TMath::Min(xmin,background->GetMinimum(v0.Data()));
   ymax = TMath::Max(ymax,background->GetMaximum(v1.Data()));
   ymin = TMath::Min(ymin,background->GetMinimum(v1.Data()));


   TH2D *hs=new TH2D("hs","",nbin,xmin,xmax,nbin,ymin,ymax);   
   TH2D *hb=new TH2D("hb","",nbin,xmin,xmax,nbin,ymin,ymax);   
   hs->SetXTitle(v0);
   hs->SetYTitle(v1);
   hb->SetXTitle(v0);
   hb->SetYTitle(v1);
   hs->SetMarkerColor(4);
   hb->SetMarkerColor(2);


   TH2F * hist = new TH2F( "MVA",    "MVA",    nbin,xmin,xmax,nbin,ymin,ymax);

   // Prepare input tree (this must be replaced by your data source)
   // in this example, there is a toy tree with signal and one with background events
   // we'll later on use only the "signal" events for the test in this example.

   Float_t MinMVA=10000, MaxMVA=-100000;
   for (Int_t ibin=1; ibin<nbin+1; ibin++){
      for (Int_t jbin=1; jbin<nbin+1; jbin++){
         var0 = hs->GetXaxis()->GetBinCenter(ibin);
         var1 = hs->GetYaxis()->GetBinCenter(jbin);
         Float_t mvaVal=reader->EvaluateMVA( "M1" ) ;
         if (MinMVA>mvaVal) MinMVA=mvaVal;
         if (MaxMVA<mvaVal) MaxMVA=mvaVal;
         hist->SetBinContent(ibin,jbin, mvaVal);
      }
   }

   // creating a fine histograms containing the error rate
   const Int_t nValBins=100;
   Double_t    sum = 0.;

   TH1F *mvaS= new TH1F("mvaS","",nValBins,MinMVA,MaxMVA);
   TH1F *mvaB= new TH1F("mvaB","",nValBins,MinMVA,MaxMVA);
   TH1F *mvaSC= new TH1F("mvaSC","",nValBins,MinMVA,MaxMVA);
   TH1F *mvaBC= new TH1F("mvaBC","",nValBins,MinMVA,MaxMVA);

   Long64_t nentries;
   nentries = TreeS->GetEntries();
//.........这里部分代码省略.........

void DrawErrors(Int_t NBINS=50, const char * filename="spin.root")
{
  TFile * infile = new TFile(filename,"READ");
  TTree * tr = (TTree*) infile->Get("asy");

  // set binning
  // -- *_div = lower limit of each bin; last one is the upper limit
  // -- *_width = bin width
  Int_t phi_bins0, eta_bins0, pt_bins0, en_bins0;
  if(gSystem->Getenv("PHI")==NULL){fprintf(stderr,"ERROR: source env vars\n"); return;};
  sscanf(gSystem->Getenv("PHI"),"%d",&phi_bins0);
  sscanf(gSystem->Getenv("ETA"),"%d",&eta_bins0);
  sscanf(gSystem->Getenv("PT"),"%d",&pt_bins0);
  sscanf(gSystem->Getenv("EN"),"%d",&en_bins0);
  const Double_t pi=3.1415;
  const Double_t phi_bins=phi_bins0; 
  const Double_t eta_bins=eta_bins0;
    const Double_t eta_low=2.6; const Double_t eta_high=4.2;
  const Double_t pt_bins=pt_bins0; 
    const Double_t pt_low=0; const Double_t pt_high=10;
  const Double_t en_bins=en_bins0; 
    const Double_t en_low=0; const Double_t en_high=100;
  const Double_t phi_low = (-1*pi)-0.1;
  const Double_t phi_high = pi+0.1;
  Double_t phi_div[phi_bins+1];
  Double_t phi_width = (phi_high - phi_low)/phi_bins;
  for(Int_t i=0; i<phi_bins; i++) phi_div[i] = phi_low + i * phi_width;
  Double_t eta_div[eta_bins+1];
  Double_t eta_width = (eta_high - eta_low)/eta_bins;
  for(Int_t i=0; i<eta_bins; i++) eta_div[i] = eta_low + i * eta_width;
  Double_t pt_div[pt_bins+1];
  Double_t pt_width = (pt_high - pt_low)/pt_bins;
  for(Int_t i=0; i<pt_bins; i++) pt_div[i] = pt_low + i * pt_width;
  Double_t en_div[en_bins+1];
  Double_t en_width = (en_high - en_low)/en_bins;
  for(Int_t i=0; i<en_bins; i++) en_div[i] = en_low + i * en_width;


  eta_div[eta_bins] = eta_high;
  pt_div[pt_bins] = pt_high;
  en_div[en_bins] = en_high;

  Float_t R3_min, R3_max;
  Float_t R3_err_min, R3_err_max;
  Float_t PB_min, PB_max;
  Float_t PB_err_min, PB_err_max;
  Float_t PY_min, PY_max;
  Float_t PY_err_min, PY_err_max;
  Float_t A_LL_min, A_LL_max;
  Float_t A_LL_err_min, A_LL_err_max;

  R3_min = tr->GetMinimum("R3");
  R3_max = tr->GetMaximum("R3");
  R3_err_min = tr->GetMinimum("R3_err");
  R3_err_max = tr->GetMaximum("R3_err");
  PB_min = tr->GetMinimum("PB");
  PB_max = tr->GetMaximum("PB");
  PB_err_min = tr->GetMinimum("PB_err");
  PB_err_max = tr->GetMaximum("PB_err");
  PY_min = tr->GetMinimum("PY");
  PY_max = tr->GetMaximum("PY");
  PY_err_min = tr->GetMinimum("PY_err");
  PY_err_max = tr->GetMaximum("PY_err");
  A_LL_min = tr->GetMinimum("A_LL");
  A_LL_max = tr->GetMaximum("A_LL");
  A_LL_err_min = tr->GetMinimum("A_LL_err");
  A_LL_err_max = tr->GetMaximum("A_LL_err");

  TH2F * R3_dist = new TH2F("R3_dist","#sigma(R_{3}) vs. R_{3}",
    NBINS,R3_min,R3_max,NBINS,R3_err_min,R3_err_max);
  TH2F * PB_dist = new TH2F("PB_dist","#sigma(P_{B}) vs. P_{B}",
    NBINS,PB_min,PB_max,NBINS,PB_err_min,PB_err_max);
  TH2F * PY_dist = new TH2F("PY_dist","#sigma(P_{Y}) vs. P_{Y}",
    NBINS,PY_min,PY_max,NBINS,PY_err_min,PY_err_max);
  TH2F * A_LL_dist[eta_bins][pt_bins][en_bins];
  TH2F * A_LL_dist_all  = new TH2F("A_LL_dist_all","#sigma(A_{LL}) vs. A_{LL} :: all bins",
    NBINS,A_LL_min,A_LL_max,NBINS,A_LL_err_min,A_LL_err_max);

  char zero_bin[128];
  strcpy(zero_bin,"phi_bin==1 && eta_bin==0 && pt_bin==0 && pt_bin==0");

  tr->Project("R3_dist","R3_err:R3",zero_bin);
  tr->Project("PB_dist","PB_err:PB",zero_bin);
  tr->Project("PY_dist","PY_err:PY",zero_bin);
  tr->Project("A_LL_dist_all","A_LL_err:A_LL");

  TFile * outfile = new TFile("error.root","RECREATE");
  R3_dist->Write();
  PB_dist->Write();
  PY_dist->Write();
  A_LL_dist_all->Write();

  char A_LL_dist_n[eta_bins][pt_bins][en_bins][256];
  char A_LL_dist_t[eta_bins][pt_bins][en_bins][256];
  char A_LL_cut[eta_bins][pt_bins][en_bins][256];
  for(Int_t g=0; g<eta_bins; g++)
  {
    for(Int_t p=0; p<pt_bins; p++)
    {
      for(Int_t e=0; e<en_bins; e++)
//.........这里部分代码省略.........

double QAnalysis::GetMaximum(string name) {
  TTree* tree = ReadTree(file);
  double max  = tree->GetMaximum(name.c_str());
  return max;
}

void nbx_step_check(const char * filename="counts.root")
{
  TFile * infile = new TFile(filename,"READ");
  TTree * tr = (TTree*) infile->Get("sca");

  Int_t NRUNS_tmp = tr->GetMaximum("i");
  const Int_t NRUNS = NRUNS_tmp;
  TH1D * h[NRUNS];
  char h_n[NRUNS][256];
  char h_t[NRUNS][256];
  char cut[NRUNS][256];

  TCanvas * c = new TCanvas("c","c",1400,1000);

  Double_t max;

  Double_t step[NRUNS][120];
  Double_t step_e[NRUNS][120];
  Double_t bx_arr[120];
  Double_t bx_arr_e[120];
  for(Int_t b=0; b<120; b++)
  {
    bx_arr[b]=b;
    bx_arr_e[b]=0;
  };
  TGraphErrors * tg[NRUNS];

  for(Int_t i=0; i<NRUNS; i++)
  {
    printf("i=%d\n",i);

    sprintf(h_n[i],"nbx_bx_%d",i+1);
    sprintf(h_t[i],"N_{bx} vs. bXing for i=%d",i+1);
    h[i] = new TH1D(h_n[i],h_n[i],120,0,120);
    sprintf(cut[i],"tot_bx*(i==%d)",i+1);
    tr->Project(h_n[i],"bx",cut[i]);

    max = h[i]->GetMaximum();
    //h[i]->Scale(1/max);




    for(Int_t b=1; b<=119; b++)
    {
      step[i][b] = h[i]->GetBinContent(b) - h[i]->GetBinContent(b+1);
      step_e[i][b] = sqrt(step[i][b]);
    };
    tg[i] = new TGraphErrors(120,bx_arr,step[i],bx_arr_e,step_e[i]);

    c->Clear();
    c->SetGrid(1,1);
    tg[i]->SetMarkerStyle(kFullCircle); 
    tg[i]->Draw("ape");

    if(i==0) c->Print("nbx_step_check.pdf(","pdf");
    else if(i+1==NRUNS) c->Print("nbx_step_check.pdf)","pdf");
    else c->Print("nbx_step_check.pdf");

  };

};

//.........这里部分代码省略.........
  int totalCand = 0;
  int totalCandInMassWindow = 0;
  int totalCandInEtaAcceptance = 0;
  int totalCandEtAbove10GeV = 0;
  int totalCandMatchedToGen = 0;

  // Loop over files
  for(UInt_t ifile=0; ifile<ntupleFileNames.size(); ifile++){

    //
    // Access samples and fill histograms
    //  
    TFile *infile = 0;
    TTree *eventTree = 0;
        
    // Data structures to store info from TTrees
    mithep::TEventInfo *info = new mithep::TEventInfo();
    mithep::TGenInfo   *gen  = new mithep::TGenInfo();
    TClonesArray *dielectronArr = new TClonesArray("mithep::TDielectron");
    
    // Read input file
    cout << "Processing " << ntupleFileNames[ifile] << "..." << endl;
    infile = new TFile(ntupleFileNames[ifile]); 
    assert(infile);
    
    // Get the TTrees
    eventTree = (TTree*)infile->Get("Events"); assert(eventTree);

    // Set branch address to structures that will store the info  
    eventTree->SetBranchAddress("Info",&info);                TBranch *infoBr       = eventTree->GetBranch("Info");

    // check whether the file is suitable for the requested run range
    UInt_t runNumMin = UInt_t(eventTree->GetMinimum("runNum"));
    UInt_t runNumMax = UInt_t(eventTree->GetMaximum("runNum"));
    std::cout << "runNumMin=" << runNumMin << ", runNumMax=" << runNumMax << "\n";
    if (!triggers.validRunRange(runNumMin,runNumMax)) {
      std::cout << "... file contains uninteresting run range\n";
      continue;
    }

    // Define other branches
    eventTree->SetBranchAddress("Dielectron",&dielectronArr); TBranch *dielectronBr = eventTree->GetBranch("Dielectron");
    TBranch *genBr = 0;
    if(sample != DATA){
      eventTree->SetBranchAddress("Gen",&gen);
      genBr = eventTree->GetBranch("Gen");
    }

    // loop over events    
    eventsInNtuple += eventTree->GetEntries();
     for(UInt_t ientry=0; ientry<eventTree->GetEntries(); ientry++) {
       //  for(UInt_t ientry=0; ientry<200000; ientry++) { // This is for faster turn-around in testing
       
      if(sample != DATA)
	genBr->GetEntry(ientry);
      
      // Check that the whole event has fired the appropriate trigger
      infoBr->GetEntry(ientry);
      
      /* Old code
      // For EPS2011 for both data and MC (starting from Summer11 production)
      // we use a special trigger for tag and probe that has second leg
      // unbiased with cuts at HLT
      ULong_t eventTriggerBit = kHLT_Ele17_CaloIdVT_CaloIsoVT_TrkIdT_TrkIsoVT_SC8_Mass30
	| kHLT_Ele32_CaloIdL_CaloIsoVL_SC17;
      // The tag trigger bit matches the "electron" of the trigger we

//.........这里部分代码省略.........
  TH2F * bz_e = new TH2F("bz_e","BBCE vs ZDCE Rate Correlation",NBINS,0,zdce_max_rate,NBINS,0,bbce_max_rate);
  TH2F * bz_w = new TH2F("bz_w","BBCW vs ZDCW Rate Correlation",NBINS,0,zdcw_max_rate,NBINS,0,bbcw_max_rate);
  TH2F * bz_x = new TH2F("bz_x","BBCX vs ZDCX Rate Correlation",NBINS,0,zdcx_max_rate,NBINS,0,bbcx_max_rate);
  TH2F * bv_e = new TH2F("bv_e","BBCE vs VPDE Rate Correlation",NBINS,0,vpde_max_rate,NBINS,0,bbce_max_rate);
  TH2F * bv_w = new TH2F("bv_w","BBCW vs VPDW Rate Correlation",NBINS,0,vpdw_max_rate,NBINS,0,bbcw_max_rate);
  TH2F * bv_x = new TH2F("bv_x","BBCX vs VPDX Rate Correlation",NBINS,0,vpdx_max_rate,NBINS,0,bbcx_max_rate);
  TH2F * vz_e = new TH2F("vz_e","VPDE vs ZDCE Rate Correlation",NBINS,0,zdce_max_rate,NBINS,0,vpde_max_rate);
  TH2F * vz_w = new TH2F("vz_w","VPDW vs ZDCW Rate Correlation",NBINS,0,zdcw_max_rate,NBINS,0,vpdw_max_rate);
  TH2F * vz_x = new TH2F("vz_x","VPDX vs ZDCX Rate Correlation",NBINS,0,zdcx_max_rate,NBINS,0,vpdx_max_rate);

  // projections
  tr->Project("bbc_ew","bbce/t:bbcw/t","bbce>0 && bbcw>0 && t>0");
  tr->Project("bbc_ex","bbce/t:bbcx/t","bbce>0 && bbcx>0 && t>0");
  tr->Project("bbc_wx","bbcw/t:bbcx/t","bbcw>0 && bbcx>0 && t>0");
  tr->Project("zdc_ew","zdce/t:zdcw/t","zdce>0 && zdcw>0 && t>0");
  tr->Project("zdc_ex","zdce/t:zdcx/t","zdce>0 && zdcx>0 && t>0");
  tr->Project("zdc_wx","zdcw/t:zdcx/t","zdcw>0 && zdcx>0 && t>0");
  tr->Project("vpd_ew","vpde/t:vpdw/t","vpde>0 && vpdw>0 && t>0");
  tr->Project("vpd_ex","vpde/t:vpdx/t","vpde>0 && vpdx>0 && t>0");
  tr->Project("vpd_wx","vpdw/t:vpdx/t","vpdw>0 && vpdx>0 && t>0");
  tr->Project("bz_e","bbce/t:zdce/t","bbce>0 && zdce>0 && t>0");
  tr->Project("bz_w","bbcw/t:zdcw/t","bbcw>0 && zdcw>0 && t>0");
  tr->Project("bz_x","bbcx/t:zdcx/t","bbcx>0 && zdcx>0 && t>0");
  tr->Project("bv_e","bbce/t:vpde/t","bbce>0 && vpde>0 && t>0");
  tr->Project("bv_w","bbcw/t:vpdw/t","bbcw>0 && vpdw>0 && t>0");
  tr->Project("bv_x","bbcx/t:vpdx/t","bbcx>0 && vpdx>0 && t>0");
  tr->Project("vz_e","vpde/t:zdce/t","vpde>0 && zdce>0 && t>0");
  tr->Project("vz_w","vpdw/t:zdcw/t","vpdw>0 && zdcw>0 && t>0");
  tr->Project("vz_x","vpdx/t:zdcx/t","vpdx>0 && zdcx>0 && t>0");
  */


  // define run total variables (total scaler counts in a single run)
  Int_t IMAX_tmp = tr->GetMaximum("i");
  const Int_t IMAX = IMAX_tmp;
  Long_t bbce_total[IMAX];
  Long_t bbcw_total[IMAX];
  Long_t bbcx_total[IMAX];
  Long_t zdce_total[IMAX];
  Long_t zdcw_total[IMAX];
  Long_t zdcx_total[IMAX];
  Long_t vpde_total[IMAX];
  Long_t vpdw_total[IMAX];
  Long_t vpdx_total[IMAX];
  Int_t time[IMAX];
  for(Int_t i=0; i<IMAX; i++)
  {
    bbce_total[i] = 0;
    bbcw_total[i] = 0;
    bbcx_total[i] = 0;
    zdce_total[i] = 0;
    zdcw_total[i] = 0;
    zdcx_total[i] = 0;
    vpde_total[i] = 0;
    vpdw_total[i] = 0;
    vpdx_total[i] = 0;
  };


  // tree loop --> fills run totals arrays
  Int_t ii;
  for(Int_t i=0; i<tr->GetEntries(); i++)
  {
    tr->GetEntry(i);
    ii = index-1; // run index starts at 1; arrays start at 0
    time[ii] = t;

void csv(TString input="tmva.csvoutput.txt", TString par1="par2", TString par2="par3", TString par3="", TString value="eventEffScaled_5") {
  std::cout << "Usage:" << std::endl
            << ".x scripts/csv.C    with default arguments" << std::endl
            << ".x scripts/csv.C(filename, par1, par2, value)" << std::endl
            << std::endl
            << "  Optional arguments:" << std::endl
            << "    filename        path to CSV file" << std::endl
            << "    par1            name of X-parameter branch" << std::endl
            << "    par2            name of Y-parameter branch (if empty, efficiency is drawn as a function of par1)" << std::endl
            << "    value           name of result (efficiency) branch" << std::endl
            << std::endl;

  TTree *tree = new TTree("data", "data");
  tree->ReadFile(input);

  gStyle->SetPalette(1);
  gStyle->SetPadRightMargin(0.14);

  TCanvas *canvas = new TCanvas("csvoutput", "CSV Output", 1200, 900);

  tree->SetMarkerStyle(kFullDotMedium);
  tree->SetMarkerColor(kRed);
  if(par2.Length() > 0) {
    //tree->Draw(Form("%s:%s", par2.Data(), par1.Data()));
    if(par3.Length() > 0)
      tree->Draw(Form("%s:%s:%s:%s", par1.Data(), par2.Data(), par3.Data(), value.Data()), "", "COLZ"); //, "", "Z");
    else
      tree->Draw(Form("%s:%s:%s", par2.Data(), par1.Data(), value.Data()), "", "COLZ"); //, "", "Z");

    TH1 *histo = tree->GetHistogram();
    if(!histo)
      return;

    histo->SetTitle(Form("%s with different classifier parameters", value.Data()));
    histo->GetXaxis()->SetTitle(Form("Classifier parameter %s", par1.Data()));
    histo->GetYaxis()->SetTitle(Form("Classifier parameter %s", par2.Data()));
    if(par3.Length() > 0)
      histo->GetZaxis()->SetTitle(Form("Classifier parameter %s", par3.Data()));
    else
      histo->GetZaxis()->SetTitle("");

    if(par3.Length() == 0) {
      float x = 0;
      float y = 0;
      float val = 0;
      double maxVal = tree->GetMaximum(value);
      double minVal = tree->GetMinimum(value);

      tree->SetBranchAddress(par1, &x);
      tree->SetBranchAddress(par2, &y);
      tree->SetBranchAddress(value, &val);
      TLatex l;
      l.SetTextSize(0.03);
    
      Long64_t nentries = tree->GetEntries();
      for(Long64_t entry=0; entry < nentries; ++entry) {
        tree->GetEntry(entry);
    
        l.SetTextColor(textColor(val, maxVal, minVal));
        l.DrawLatex(x, y, Form("%.3f", val*100));
      }
    }
  }
  else {
    tree->Draw(Form("%s:%s", value.Data(), par1.Data()));
    TH1 *histo = tree->GetHistogram();
    if(!histo) 
      return;
    histo->SetTitle(Form("%s with different classifier parameters", value.Data()));
    histo->GetXaxis()->SetTitle(Form("Classifier parameter %s", par1.Data()));
    histo->GetYaxis()->SetTitle(value);
  }
}

void mk_tree(const char * acc_file="datfiles/acc.dat")
{
  // read acc file into tree
  TTree * acc = new TTree("acc","counts tree from acc.dat");
  char cols[2048];
  char bbc_cols[256];
  char zdc_cols[256];
  char vpd_cols[256];
  for(Int_t i=0; i<=7; i++)
  {
    if(i==0) 
    {
      sprintf(bbc_cols,"bbc_%d/D",i);
      sprintf(zdc_cols,"zdc_%d/D",i);
      sprintf(vpd_cols,"vpd_%d/D",i);
    }
    else
    {
      sprintf(bbc_cols,"%s:bbc_%d/D",bbc_cols,i);
      sprintf(zdc_cols,"%s:zdc_%d/D",zdc_cols,i);
      sprintf(vpd_cols,"%s:vpd_%d/D",vpd_cols,i);
    };
  };
  sprintf(cols,"i/I:runnum/I:fi/I:fill/I:t/D:bx/I:%s:%s:%s:tot_bx/D:blue/I:yell/I",bbc_cols,zdc_cols,vpd_cols);
  printf("%s\n",cols);
  acc->ReadFile(acc_file,cols);

  acc->Print();

  Int_t IMAX_tmp = acc->GetMaximum("i");
  const Int_t IMAX = IMAX_tmp;

  // set branch addresses to read through acc tree
  Int_t index,runnum,fill_index,fill,bx;
  Double_t bbc[8];
  Double_t zdc[8];
  Double_t vpd[8];
  Double_t time;
  Double_t tot_bx;
  Int_t blue,yell;
  acc->SetBranchAddress("i",&index);
  acc->SetBranchAddress("runnum",&runnum);
  acc->SetBranchAddress("fi",&fill_index);
  acc->SetBranchAddress("fill",&fill);
  acc->SetBranchAddress("t",&time);
  acc->SetBranchAddress("bx",&bx);
  char str[16];
  for(Int_t i=0; i<8; i++) { sprintf(str,"bbc_%d",i); acc->SetBranchAddress(str,&bbc[i]); };
  for(Int_t i=0; i<8; i++) { sprintf(str,"zdc_%d",i); acc->SetBranchAddress(str,&zdc[i]); };
  for(Int_t i=0; i<8; i++) { sprintf(str,"vpd_%d",i); acc->SetBranchAddress(str,&vpd[i]); };
  acc->SetBranchAddress("tot_bx",&tot_bx);
  acc->SetBranchAddress("blue",&blue);
  acc->SetBranchAddress("yell",&yell);


  // build arrays for restructuring; arrays are needed so that
  // we can implement bXing shift corrections
  Double_t bbce_arr[IMAX][120];
  Double_t bbcw_arr[IMAX][120];
  Double_t bbcx_arr[IMAX][120];
  Double_t zdce_arr[IMAX][120];
  Double_t zdcw_arr[IMAX][120];
  Double_t zdcx_arr[IMAX][120];
  Double_t vpde_arr[IMAX][120];
  Double_t vpdw_arr[IMAX][120];
  Double_t vpdx_arr[IMAX][120];
  Int_t runnum_arr[IMAX];
  Int_t fi_arr[IMAX];
  Int_t fill_arr[IMAX];
  Double_t time_arr[IMAX];
  Double_t tot_bx_arr[IMAX][120];
  Int_t blue_arr[IMAX][120];
  Int_t yell_arr[IMAX][120];
  Bool_t kicked_arr[IMAX][120];


  // restructure tree into one suitable for analysis
  TTree * sca = new TTree("sca","restructured tree");
  Double_t bbce,bbcw,bbcx; // e=east, w=west, x=coincidence
  Double_t zdce,zdcw,zdcx;
  Double_t vpde,vpdw,vpdx;
  Bool_t okEntry,kicked;
  sca->Branch("i",&index,"i/I");
  sca->Branch("runnum",&runnum,"runnum/I");
  sca->Branch("fi",&fill_index,"fi/I");
  sca->Branch("fill",&fill,"fill/I");
  sca->Branch("t",&time,"t/D");
  sca->Branch("bx",&bx,"bx/I");
  sca->Branch("bbce",&bbce,"bbce/D");
  sca->Branch("bbcw",&bbcw,"bbcw/D");
  sca->Branch("bbcx",&bbcx,"bbcx/D");
  sca->Branch("zdce",&zdce,"zdce/D");
  sca->Branch("zdcw",&zdcw,"zdcw/D");
  sca->Branch("zdcx",&zdcx,"zdcx/D");
  sca->Branch("vpde",&vpde,"vpde/D");
  sca->Branch("vpdw",&vpdw,"vpdw/D");
  sca->Branch("vpdx",&vpdx,"vpdx/D");
  sca->Branch("tot_bx",&tot_bx,"tot_bx/D");
  sca->Branch("blue",&blue,"blue/I");
  sca->Branch("yell",&yell,"yell/I");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      sprintf(vpd_br[i],"vpd_%d",i);
      acc->SetBranchAddress(vpd_br[i],&(vpd[i]));
    //};
  };

  // ----------------------------------------------

  // total counts vs. scaler bit bar charts
  TH1F * ntot_vs_bits_bbc = new TH1F();
  TH1F * ntot_vs_bits_zdc = new TH1F();
  TH1F * ntot_vs_bits_vpd = new TH1F();


  // scaler bit combination names
  char comb[8][16];
  strcpy(comb[0],"none");
  strcpy(comb[1],"e");
  strcpy(comb[2],"w");
  strcpy(comb[3],"w+e");
  strcpy(comb[4],"x");
  strcpy(comb[5],"x+e");
  strcpy(comb[6],"x+w");
  strcpy(comb[7],"x+w+e");

  // fill bar charts
  for(Int_t i=0; i<acc->GetEntries(); i++)
  {
    acc->GetEntry(i);
    for(Int_t j=0; j<8; j++)
    {
      ntot_vs_bits_bbc->Fill(comb[j],bbc[j]);
      ntot_vs_bits_zdc->Fill(comb[j],zdc[j]);
      //if(j<4) ntot_vs_bits_vpd->Fill(comb[j],vpd[j]);
      ntot_vs_bits_vpd->Fill(comb[j],vpd[j]);
    };
  };

  ntot_vs_bits_bbc->SetStats(0);
  ntot_vs_bits_bbc->SetTitle("total bbc counts vs. scaler bits");
  ntot_vs_bits_bbc->SetBarWidth(0.4);
  ntot_vs_bits_bbc->SetBarOffset(0.55);
  ntot_vs_bits_bbc->SetFillColor(50);
  TCanvas * c_bbc_bits = new TCanvas("c_bbc_bits","c_bbc_bits",700,500);
  c_bbc_bits->SetGrid(0,1);
  c_bbc_bits->SetLogy();
  ntot_vs_bits_bbc->Draw("bar2");

  ntot_vs_bits_zdc->SetStats(0);
  ntot_vs_bits_zdc->SetTitle("total zdc counts vs. scaler bits");
  ntot_vs_bits_zdc->SetBarWidth(0.4);
  ntot_vs_bits_zdc->SetBarOffset(0.55);
  ntot_vs_bits_zdc->SetFillColor(50);
  TCanvas * c_zdc_bits = new TCanvas("c_zdc_bits","c_zdc_bits",700,500);
  c_zdc_bits->SetGrid(0,1);
  c_zdc_bits->SetLogy();
  ntot_vs_bits_zdc->Draw("bar2");

  ntot_vs_bits_vpd->SetStats(0);
  ntot_vs_bits_vpd->SetTitle("total vpd counts vs. scaler bits");
  ntot_vs_bits_vpd->SetBarWidth(0.4);
  ntot_vs_bits_vpd->SetBarOffset(0.55);
  ntot_vs_bits_vpd->SetFillColor(50);
  TCanvas * c_vpd_bits = new TCanvas("c_vpd_bits","c_vpd_bits",700,500);
  c_vpd_bits->SetGrid(0,1);
  c_vpd_bits->SetLogy();
  ntot_vs_bits_vpd->Draw("bar2");

  char outdir[32];
  char bbc_outfile[64];
  char zdc_outfile[64];
  char vpd_outfile[64];
  strcpy(outdir,"bit_combos");
  sprintf(bbc_outfile,"%s/bbc_bit_combos.png",outdir);
  sprintf(zdc_outfile,"%s/zdc_bit_combos.png",outdir);
  sprintf(vpd_outfile,"%s/vpd_bit_combos.png",outdir);
  c_bbc_bits->Print(bbc_outfile,"png");
  c_zdc_bits->Print(zdc_outfile,"png");
  c_vpd_bits->Print(vpd_outfile,"png");
  printf("%s created\n",bbc_outfile);
  printf("%s created\n",zdc_outfile);
  printf("%s created\n",vpd_outfile);


  // ----------------------------------------------


  // get maximum number of runs  = IMAX
  Int_t IMAX_tmp = acc->GetMaximum("i");
  const Int_t IMAX = IMAX_tmp;


  // compute no. bXings with possible interaction in a given run
  Int_t nbx[IMAX]; 
  for(Int_t i=0; i<IMAX; i++) nbx[i]=0;
  for(Int_t i=0; i<acc->GetEntries(); i++)
  {
    acc->GetEntry(i);
    if(blue*yell != 0) nbx[index-1]++;
  };
};

void makePlotWithSelection(TTree* tr, TString varToPlot, TString canvName, TString xAxisName, TCut cutLoose, TCut cutAdditional)
{

  TTree* trCutted = tr->CopyTree(cutLoose);
  TH1D* histLoose = new TH1D("hist",varToPlot,50,trCutted->GetMinimum(varToPlot)-0.1*fabs(trCutted->GetMinimum(varToPlot)),trCutted->GetMaximum(varToPlot)+0.1*fabs(trCutted->GetMaximum(varToPlot)));
  trCutted->Draw(varToPlot+TString(">>hist"),"1","goff");
  TCanvas* canv = new TCanvas(canvName,canvName);
  //TAxis* xAxis = histLoose->GetXaxis();
  histLoose->GetXaxis()->SetTitle(xAxisName);
  histLoose->GetYaxis()->SetRangeUser(0,histLoose->GetMaximumStored());
  histLoose->SetFillStyle(3004);
  histLoose->SetFillColor(2);
  histLoose->Draw();
  trCutted->SetFillStyle(1001);
  trCutted->SetFillColor(3);
  trCutted->Draw(varToPlot,cutAdditional,"same");
  canv->SaveAs(canvName+".png");
  delete histLoose;
  histLoose=0;
  trCutted=0;
  //xAxis=0;
}