C++ TH2D::Draw方法代码示例

void SimplePlot2D(){

  //string filename = "data/scan_CH1-64_unmasked.root";
  //string filename = "data/scan_CH1-50_masked.root";
  string filename = "test.root";
  string varXname = "VMM #";
  string varYname = "CH #";
  
  // delay count stuff
  int CH = 21;
  double delays[5];
  double count_tot[5];
  double count_right[5];
  for(int i = 0; i < 5; i++){
    delays[i] = double(i)*5.;
    count_tot[i] = 0.;
    count_right[i] = 0.;
  }
  
      

  ///////////////////////////////////////////////////////
  
  TChain* tree = new TChain("MMFE8","MMFE8");

  tree->AddFile(filename.c_str());

  MMFE8Base* base = new MMFE8Base(tree);

  int N = tree->GetEntries();

  TH2D* hist = new TH2D("hist","hist", 8, 0.5, 8.5, 64, 0.5,64.5);
  TH2D* histN = (TH2D*) hist->Clone("norm");
  TH2D* histchch = new TH2D("histchch","histchch", 64, 0.5, 64.5, 64, 0.5,64.5);
  TH2D* histDelay = new TH2D("histN","histN", 31, 9.5, 40.5, 5,-0.5, 4.5);	
  TH2D* histDelayD = new TH2D("histD","histD", 31, 9.5, 40.5, 5,-0.5, 4.5);

  for(int i = 0; i < N; i++){
    base->GetEntry(i);
    
    if(base->CHpulse == CH){
      //count_tot[base->Delay] += 1.;
      count_tot[(base->TPDAC-80)/20] += 1.;
      if(base->CHpulse == base->CHword)
	//count_right[base->Delay] += base->TDO;
      	count_right[(base->TPDAC-80)/20] += base->PDO;
    }

    //histDelayD->Fill(base->CHpulse,base->Delay);
    histDelayD->Fill(base->CHpulse,(base->TPDAC-80)/20);
    if(base->CHpulse == base->CHword)
      //histDelay->Fill(base->CHpulse,base->Delay,base->TDO);
      histDelay->Fill(base->CHpulse,(base->TPDAC-80)/20,base->PDO);

    if((base->CHpulse != base->CHword || true) &&
       base->VMM == 6)
      histchch->Fill(base->CHpulse,base->CHword);

    if(base->CHpulse != base->CHword)
      continue;

    hist->Fill(base->VMM,base->CHpulse,base->PDO);
    histN->Fill(base->VMM,base->CHpulse);
  
  }

  for(int x = 0; x < 8; x++){
    for(int y = 0; y < 64; y++){
      double v = hist->GetBinContent(x+1,y+1);
      double N = histN->GetBinContent(x+1,y+1);
      hist->SetBinContent(x+1,y+1,v/max(int(N),1));
    }
  }
  
  TLatex l;
  //l.NDC();

  TCanvas* can = new TCanvas("can","can",600,500);
  can->SetTopMargin(0.05);
  can->SetLeftMargin(0.12);
  can->SetRightMargin(0.15);
  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);

  can->Draw();
  can->SetGridx();
  can->SetGridy();
  
  can->cd();

  hist->Draw("COLZ");
  
  hist->GetXaxis()->SetTitle(varXname.c_str());
  hist->GetXaxis()->CenterTitle();
  hist->GetYaxis()->SetTitle(varYname.c_str());
  hist->GetYaxis()->CenterTitle();
  hist->GetYaxis()->SetTitleOffset(1.4);
  hist->GetYaxis()->CenterTitle();
  //hist->GetYaxis()->SetRangeUser(0.,hist->GetMaximum()*1.1) ;

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

// infilename - root file with relevant histograms
// system - PP,APAP,PP
// status - Pass,Fail
// rWrite - 0-no,1-png,2-eps
// rPerformance - 0-no,1-yes (ALICE logo etc.)
// bin: 0 - all, 1- 0:5, 2- 5:10, etc
void drawPID(const char* infilename, const char* system, const char* status, Int_t rWrite, Int_t rPerformance, Int_t bin)
{


    TFile *f = new TFile(infilename, "read");

    // TPC dEdx
    TH2D* TPCdEdx =(TH2D*)f->Get(Form("TPCdEdxcut%s1%stpcM%i",status, system,0));

    if (!bin) {
        int minMultBin = 0;
        int maxMultBin = 6; // 8
    }
    else {
        int minMultBin = bin-1;
        int maxMultBin = bin; // 8
    }

    double EvMultall = 0;

    for (int i = minMultBin; i < maxMultBin; i++) {

        TH2D* TPCdEdxN =(TH2D*)f->Get(Form("TPCdEdxcut%s1%stpcM%i",status,system,i));
        TPCdEdx->Add(TPCdEdxN);

        cout << i << " " << TPCdEdxN->GetEntries() << endl;

        //delete hEvMult;
    }

    TCanvas *c2 = new TCanvas("TPC dEdx", "TPC dEdx");
    c2->SetGridx();
    c2->SetGridy();
    c2->SetFillColor(10);
    c2->SetRightMargin(1.9);
    c2->SetLogz();

    TPCdEdx->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
    TPCdEdx->GetXaxis()->SetRangeUser(0.0,6.0);
    TPCdEdx->GetYaxis()->SetTitle("dE/dx");
    TPCdEdx->GetZaxis()->SetLabelSize(0.03);
    TPCdEdx->Draw("colz");

// 	double a1 = -3000.0; double b1 =  1280.0;
// 	double a2 = -312.5;  double b2 =  312.5;
// 	double a3 = -200.0;  double b3 =  240.0;

//    TF1 *fa1 = new TF1("fa1","-1800*x+940",0.3,0.4);
//   fa1->Draw("same");
//   TF1 *fa2 = new TF1("fa2","-500.0*x+420.0",0.4,0.6);
//   fa2->Draw("same");
// TF1 *fa3 = new TF1("fa3","-216.7*x+250.0",0.6,0.9);
//    fa3->Draw("same");
// TF1 *fa4 = new TF1("fa4","-566.7*x+570.0",0.6,0.75);
//    fa4->Draw("same");
// TF1 *fa5 = new TF1("fa5","-2076.92*x+1476.15",0.47,0.6);
//    fa5->Draw("same");
//   cout<<TPCdEdx->GetNbinsX()<<endl;
//   cout<<TPCdEdx->GetNbinsY()<<endl;
//   for (int ii=0;ii<TPCdEdx->GetNbinsX();ii++){

//     for (int jj=0;jj<TPCdEdx->GetNbinsY();jj++){

//       cout<<"binX: "<<ii<<endl;
//       cout<<"binY: "<<jj<<endl;
//       cout<<"val: "<<TPCdEdx->GetBinContent(ii,jj)<<endl;
//     }
//   }

// 	TH1D *py = TPCdEdx->ProjectionY("py", 230, 232); // where firstYbin = 0 and lastYbin = 9
// 	TCanvas *c22 = new TCanvas("TPC2", "TPC2");
// 	py->Draw();

    postprocess(c2,Form("TPCdEdx%s",status),rWrite,rPerformance,system);

    // TPC Nsigma
    TH2D* TPCNsigma =(TH2D*)f->Get(Form("TPCNSigmacut%s1%stpcM%i",status,system,0));

    // int minMultBin = 0;
    // int maxMultBin = 2;
    double EvMultall = 0;

    for(int i = minMultBin; i<maxMultBin; i++) {
        //all
        TH2D* TPCNsigmaN =(TH2D*)f->Get(Form("TPCNSigmacut%s1%stpcM%i",status,system,i));
        TPCNsigma->Add(TPCNsigmaN);
        //delete hEvMult;
    }

    TCanvas *c3 = new TCanvas("TPC Nsigma", "TPC Nsigma");
    c3->SetGridx();
    c3->SetGridy();
    c3->SetFillColor(10);
    c3->SetRightMargin(1.7);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      if (factor / background < 10) continue;

      // Get the total number of events in the gaussian 1 sigma
      Int_t sigRangeLow = hist->FindBin(mean - sigma);
      Int_t sigRangeHigh = hist->FindBin(mean + sigma);

      double sigRangeTotal = 0;
      for (Int_t bin = sigRangeLow; bin <= sigRangeHigh; bin++)
        sigRangeTotal += hist->GetBinContent(bin);

      // 2 * 1 sigma integral shoudl give ~ area under gaussian
      sigma /= sqrt(2 * sigRangeTotal);

      ptsX.push_back(hist->GetYaxis()->GetBinCenter(row));
      ptsY.push_back(mean);
      ptsErr.push_back(sigma);
    }

    if (ptsX.size() < 3) continue;

    std::vector<double> yvals = ptsY;
    std::sort(yvals.begin(), yvals.end());
    const double median = yvals[yvals.size()/2];
    double avgDeviation = 0;
    for (unsigned int i = 0; i < yvals.size(); i++)
      avgDeviation += fabs(yvals[i] - median);
    avgDeviation /= (double)yvals.size();

    std::vector<double> ptsXGood;
    std::vector<double> ptsYGood;
    std::vector<double> ptsErrGood;

    for (unsigned int i = 0; i < ptsX.size(); i++)
    {
      if (fabs(ptsY[i] - median) > 1.5*avgDeviation) continue;
      ptsXGood.push_back(ptsX[i]);
      ptsYGood.push_back(ptsY[i]);
      ptsErrGood.push_back(ptsErr[i]);
    }

    if (ptsXGood.size() < 3) continue;

    TGraphErrors* graph = new TGraphErrors(ptsXGood.size(),
                                           &(ptsXGood.at(0)),
                                           &(ptsYGood.at(0)), 0,
                                           &(ptsErrGood.at(0)));

    TF1* fitFunc = new TF1("f1", "1 ++ x");
    TF1* result = 0;

    graph->Fit(fitFunc, "Q0E").Get();
    result = graph->GetFunction(fitFunc->GetName());

    // Weight the angle by the slope uncertainty and the inverse of the chi2 normalized
    double weight = result->GetParError(1);
    const double chi2 = result->GetChisquare() / (double)result->GetNDF();
    fitChi2 += chi2;
    weight *= chi2;
    if (weight > 10 * DBL_MIN) weight = 1.0 / weight;
    else weight = 1.0;

    if (axis)
    {
      rotation -= weight * atan(result->GetParameter(1));
      offsetX = result->GetParameter(0);
    }
    else
    {
      rotation += weight * atan(result->GetParameter(1));
      offsetY = result->GetParameter(0);
    }

    angleWeights += weight;

    if (display)
    {
      TCanvas* can = new TCanvas("ResidualAlignment", "Residual Alignment", 900, 600);
      can->Divide(2);
      can->cd(1);
      hist->Draw("COLZ");
      can->cd(2);
      result->SetLineColor(46);
      result->SetLineWidth(2);
      graph->Draw("ap");
      result->Draw("SAME");
      can->Update();
      can->WaitPrimitive();
    }

    delete fitFunc;
    delete graph;
  }

  if (angleWeights > 10 * DBL_MIN)
    rotation /= angleWeights;
  std::cout << "relaxation: " << relaxation << std::endl;
  rotation *= relaxation;
  offsetX *= relaxation;
  offsetY *= relaxation;
}

//.........这里部分代码省略.........
    }
    delete hsignal0[i];
    delete hsignal1[i];
    delete hsignal[i];
    delete hbackground0[i];
    delete hbackground1[i];
    delete hbackground[i];
  }

  for(int nbin=1;nbin<4;nbin++)
  {
    for(int i=0;i<ntrgbins1;i++)
    {
      rn[i]=Vn[i][nbin];
      rn_err[i]=VnError[i][nbin];
    }

    gr_v2[jj][nbin] = new TGraphErrors(ntrgbins1,eta,rn,eta_err,rn_err);
    gr_v2[jj][nbin]->SetMarkerColor(color[nbin]);
  }
}

  TCanvas* c = new TCanvas("c","c",900,500);
  makeMultiPanelCanvas(c,4,2,0.01,0.0,0.2,0.2,0.02);
  TH2D* htmp = new TH2D("htmp",";#eta^{a};r_{2}(#eta^{a},#eta^{b})",100,-2.4,2.4,100,0.0,0.2);
  fixedFontHist(htmp,1.8,2.2);
  htmp->GetXaxis()->CenterTitle();
  htmp->GetYaxis()->CenterTitle();
  htmp->GetYaxis()->SetTitleSize(htmp->GetYaxis()->GetTitleSize()*1.2);
  htmp->GetXaxis()->SetTitleSize(htmp->GetXaxis()->GetTitleSize()*1.);
  for(int jj=0;jj<8;jj++)
  {
    c->cd(jj+1);
    htmp->Draw();
    gr_v2[jj][1]->SetMarkerStyle(21);
    gr_v2[jj][1]->Draw("Psame");
  }

return;
  TString histtitle[8] = {"0-0.2%","0-5%","5-10%","10-20%","20-30%","30-40%","40-50%","50-60%"};
  TLine* ll = new TLine(0,1.,5.0,1.);
  TCanvas* c = new TCanvas("c","c",900,500);
  makeMultiPanelCanvas(c,4,2,0.01,0.0,0.2,0.2,0.02);
  TH2D* htmp = new TH2D("htmp",";#eta^{a};r_{2}(#eta^{a},#eta^{b})",100,-0.1,2.23,100,0.801,1.06);
  fixedFontHist(htmp,1.8,2.2);
  htmp->GetXaxis()->CenterTitle();
  htmp->GetYaxis()->CenterTitle();
  htmp->GetYaxis()->SetTitleSize(htmp->GetYaxis()->GetTitleSize()*1.2);
  htmp->GetXaxis()->SetTitleSize(htmp->GetXaxis()->GetTitleSize()*1.);
  for(int jj=0;jj<8;jj++)
  {
    c->cd(jj+1);
    htmp->Draw();
    ll->Draw("Lsame");
    gr[jj][1]->SetMarkerStyle(21);
    gr[jj][1]->Draw("Psame");
    fit[jj][1]->Draw("Lsame");

    gr[jj+8][1]->SetMarkerStyle(25);
//    fit[jj+8][1]->SetLineStyle(5);
    gr[jj+8][1]->Draw("Psame");
//    fit[jj+8][1]->Draw("Lsame");
  }

  TLatex* latex2 = new TLatex();
  latex2->SetNDC();

int main(int argc, char * argv[])
{
    int ChooseCamera;
    cout << "Set the camera you want to use: \n 0 -> Integrated one \n 1 -> USB one " << endl;
    cin >> ChooseCamera;
    VideoCapture cam(ChooseCamera); // with 1 -> choose the webcam if it's attached, otherwise choose the integrated one  
    
    if(!cam.isOpened()){  // Controllo di apertura della webcam
        cout << "Impossibile aprire lo stream della webcam. Il programma verrà terminato." << endl;
        return -1;
	}   
        
    //~ double MyWidth = 640; // Limiti 1280x720 copiati da Ave. Sono quelli propri della videocamera come mostrati da Cheese.
    //~ double MyHeight = 480;
    double MyWidth = 1280; // Limiti 1280x720 copiati da Ave. Sono quelli propri della videocamera come mostrati da Cheese.
    double MyHeight = 720;
    cam.set(CV_CAP_PROP_FRAME_WIDTH, MyWidth);  //Set frame width and height
    cam.set(CV_CAP_PROP_FRAME_HEIGHT, MyHeight);
    //~ system("v4l2-ctl --set-fmt-video=width=1280,height=720,pixelformat=1 -d /dev/video1"); //Set frame width and height through a system call
	cam.set(CV_CAP_PROP_FPS, 30); //Set the webcam frame rate
	
	// NB: per poter modificare i valori della frame Width and Height è necessario compilare le librerie OPENCV con
	// le librerie [ v4l2ucp, v4l-utils and libv4l-dev ] già installate nel computer 
	// (vedi: http://stackoverflow.com/questions/16287488/runtime-opencv-highgui-error-highgui-error-v4l-v4l2-vidioc-s-crop-opencv-c)
	
	
	if (ChooseCamera==1) {
		//~ system("v4l2-ctl --set-ctrl exposure_auto=1 -d /dev/video1"); 			//Alternative method to change v4l2 controls
		//~ system("v4l2-ctl --set-ctrl exposure_absolute=5 -d /dev/video1");
		system("v4l2-ctl -c exposure_auto=1 -d /dev/video1");				//system call to adjust the balance variables, sono necessarie le librerie v4l2
		system("v4l2-ctl -c exposure_absolute=5 -d /dev/video1");			//Le chiamate e le impostazioni le ho copiate da Ave, sperando che siano quelle giust
		// in general v4l2-ctl -c command -> comando per cambiare le immpostazioni, lista delle impostazioni possibili le leggi 
		// da bash dando $v4l2-ctl --list-ctrls 
		// documentazione completa delle v4l2 la trovi @ http://ivtvdriver.org/index.php/V4l2-ctl
		cout << "Sys call for camera #" << ChooseCamera << " done!\n";
		//system("say 'settings done baby'");
	}

	if (ChooseCamera==0) {
		//~ //Put here commands for set the integrated camera commands. Generally useless
		//~ system("v4l2-ctl -c exposure_auto=1 -d /dev/video0");				//system call to adjust the balance variables, sono necessarie le librerie v4l2
		//~ system("v4l2-ctl -c exposure_absolute=5 -d /dev/video0");			//Le chiamate e le impostazioni le ho copiate da Ave, sperando che siano quelle giust
		cout << "Sys call for camera #" << ChooseCamera << " done!\n";
	}


	//Controllo dimensionamento immagine acquisita:
	double RealWidth = cam.get(CV_CAP_PROP_FRAME_WIDTH); 				//This function get the Width (in pixels) of the output.
	double RealHeight = cam.get(CV_CAP_PROP_FRAME_HEIGHT);				//This function get the Height (in pixels) of the output.
	cout << " Setted Width: " << MyWidth << "		|| Real Width: " << RealWidth  << endl;
	cout << " Setted Height: " << MyHeight << "		|| Real Height: " << RealHeight  << endl;
	cout << "\n";
	
	
	//Creazione Root object used for analysis
	TApplication myapp("myapp", &argc,argv);
	
	//ROOT Canvas for the current image
	//~ TCanvas * C = new TCanvas ("C","",710,400);
	//~ TH2D * CurrentHisto = new TH2D ("h2d", " ", (int)MyWidth, -(MyWidth/2)+1 , (MyWidth/2), (int)MyHeight, -(MyHeight/2)+1 , (MyHeight/2) );
	//Canvas for the total image 
	TCanvas * C2 = new TCanvas ("C2","",710,400);
	TH2D * TotalHisto = new TH2D ("TotalHisto", " ", (int)MyWidth, -(MyWidth/2)+1 , (MyWidth/2), (int)MyHeight, -(MyHeight/2)+1 , (MyHeight/2) );
	//~ gStyle->SetPalette(1);
	TotalHisto->Draw("CONT");
	
	//~ double RootTempMatrix [int(MyWidth)][int(MyHeight)]; //declaration of a two dimansional array
	//~ int meanX;
	//~ int meanY;
	//~ double VarX;
	//~ double VarY;
	//~ double IntIntegral; 
	double alfa= 0.79; //IntensityCalibration


    namedWindow("Live_Video",WINDOW_NORMAL); 			//Crea una finestra di cui posso modificare le dimensioni (openCV method)
    resizeWindow("Live_Video", RealWidth, RealHeight); 	//Imposta la dimensione iniziale uguale alla dimensione del frame
    int k=0;
    float DisplayHeight = RealHeight / 2 ;
	float DisplayWidth = RealWidth / 2 ;
    Mat ImageBGR(Size(RealWidth, RealHeight), CV_16UC3);
    
    Vec3b intensity;
	//Definisco due array sui quali metto l'immagine istantanea e l'immagine totale
    double * MatIn = new double [int( MyWidth * MyHeight)];
    double * MatSum = new double [int( MyWidth * MyHeight)];
    
    for (int m=0; m < (MyWidth * MyHeight); m++){
		MatSum[m]=0.0;
		}
	
	for(int i=0; i< int(MyWidth); i++){			//Piazzo dentro i valori
		for(int j=0; j< int(MyHeight); j++){	
			TotalHisto->SetBinContent((MyWidth-i+1),(MyHeight-j+1),0);
		}
	}
        
	int IntensityGreen;
	double temp = 0;
    
//.........这里部分代码省略.........

void DrawKinvarPlot(TString filename="kinvarset/FMSOR_Pt_vs_run.root",
                    TString classname="pi0") {
  TFile * infile = new TFile(filename.Data(),"READ");

  char trig[32];
  char typ[32];
  char tmp[128];
  TString tmpstr = filename;
  tmpstr.ReplaceAll("_"," ");
  TRegexp re("^.*\/");
  tmpstr(re) = "";
  cout << tmpstr << endl;
  sscanf(tmpstr.Data(),"%s %s",trig,typ);
  printf("%s %s\n",trig,typ);

  TString hname = "h2_"+classname;
  TH2D * h = (TH2D*)infile->Get(hname.Data());
  TString gname = "g_"+classname;
  TGraphErrors * g = (TGraphErrors*)infile->Get(gname.Data());
  TString gtname = "gt_"+classname;
  TGraphErrors * gt;
  if(!strcmp(typ,"Pt")) gt = (TGraphErrors*)infile->Get(gtname.Data());
  
  char htitle[256];
  char hnewtitle[512];
  strcpy(htitle,h->GetTitle());
  sprintf(hnewtitle,"%s -- %s triggers",htitle,trig);
  h->SetTitle(hnewtitle);

  h->SetMinimum(0.001);

  if(!strcmp(typ,"Pt")) {
    gt->SetLineWidth(2);
    gt->SetLineColor(kBlack);
  };

  // temporary hack to cut out fluctuations of pt_thresh vs. run to 0
  Double_t xx,yy;
  Int_t gtn = gt->GetN();
  for(int ii=0; ii<gtn; ii++) {
    gt->GetPoint(ii,xx,yy);
    printf("xx=%f yy=%f\n",xx,yy);
    if(yy<0.001) gt->RemovePoint(ii);
  };
   
  // drawing range
  Float_t xmin = h->GetXaxis()->GetXmin();
  Float_t xmax = h->GetXaxis()->GetXmax();
  xmax = 515;
  if(xmax<gtn) { 
    fprintf(stderr,"\nWARNING: xmax < gt->GetN(); graphs will be cut off!\n\n");
  };
  
  TCanvas * c = new TCanvas("c","c",1500,700);
  c->Divide(1,2);
  gStyle->SetOptStat(0);
  for(int x=1; x<=2; x++) c->GetPad(x)->SetGrid(1,1);
  c->cd(1);
  c->GetPad(1)->SetLogz();
  h->GetXaxis()->SetRangeUser(xmin,xmax);
  h->Draw("colz");
  if(!strcmp(typ,"Pt")) gt->Draw("LX");
  c->cd(2);
  g->SetFillColor(kGray);
  printf("%f %f\n",xmin,xmax);
  g->GetXaxis()->SetLimits(xmin,xmax);
  g->GetXaxis()->SetRangeUser(xmin,xmax);
  /*
  g->Draw("A3");
  g->Draw("PLX");
  */
  g->Draw("APLX");
  
  TString outname = filename.ReplaceAll(".root"," "+classname);
  outname = outname+".png";
  c->Print(outname.Data(),"png");
};

void RayTrace(AOpticsManager* manager, TCanvas* can3D)
{
  const int kNdeg = 8;
  TH2D* h2[kNdeg];
  TGraph* graph = new TGraph();
  TCanvas* can = new TCanvas("can", "can", 900, 900);
  TCanvas* can2= new TCanvas("can2", "can2", 900, 900);
  can->Divide(3, 3, 1e-10, 1e-10);

  TH2D* hMirror = new TH2D("hMirror", ";X (mm);Y (mm)", 1000, -7, 7, 1000, -7, 7);

  for(int i = 0; i < kNdeg; i++){

    double deg = i*0.5;
    TGeoTranslation raytr("raytr", -2*kF*TMath::Sin(deg*TMath::DegToRad()), 0, 2*kF*TMath::Cos(deg*TMath::DegToRad()));
    TVector3 dir;
    dir.SetMagThetaPhi(1, TMath::Pi() - deg*TMath::DegToRad(), 0);
    double lambda = 400*nm; // dummy
    ARayArray* array = ARayShooter::Square(lambda, 14*m, 401, 0, &raytr, &dir);

    manager->TraceNonSequential(*array);

    h2[i] = new TH2D("", Form("#it{#theta} = %3.1f#circ;x (mm); y (mm)", deg), 200, -40, 100, 200, -70, 70);
    TH2D tmp("", "", 100, -1e5, 1e5, 100, -1e5, 1e5);

    TObjArray* focused = array->GetFocused();

    for(Int_t j = 0; j <= focused->GetLast(); j++){
      ARay* ray = (ARay*)(*focused)[j];
      Double_t p[4];
      ray->GetLastPoint(p);
      tmp.Fill(p[0], p[1]);

      if (i == 0) {
        int n = ray->FindNodeNumberStartWith("mirror");
        const double* pn = ray->GetPoint(n);
        hMirror->Fill(pn[0]/m, pn[1]/m);
      } // if

      if (i == kNdeg - 1 && gRandom->Uniform() < 0.001) {
        TPolyLine3D* pol = ray->MakePolyLine3D();
        pol->SetLineColor(2);
        can3D->cd();
        pol->Draw();
      } // if
    } // j

    double meanx = tmp.GetMean();

    for(Int_t j = 0; j <= focused->GetLast(); j++){
      ARay* ray = (ARay*)(*focused)[j];
      Double_t p[4];
      ray->GetLastPoint(p);
      h2[i]->Fill((p[0] - meanx)/mm, p[1]/mm);
    } // j

    can->cd(i + 1);
    h2[i]->Draw("colz");

    if(i == 0){
      can2->cd();
      hMirror->Draw("colz");
    } // i

    delete array;
  } // i
}

int main(int argc, char** argv){
	std::cout << "Starting application" << std::endl;
	TApplication* app = new TApplication("Analyse",&argc,argv);

	//Set plot style for own root design
	PlotStyle* s = new PlotStyle();

	CanvasManager* cManager = new CanvasManager("combinedPlots");

	//#######################
	//# Creating Trees and Analyse objects
	//#######################

	std::vector<AnalysisWrapper*> wrapperVect;

	wrapperVect.push_back(new AnalysisWrapper("deltaPhiGun/DeltaPhiGunPt100dPhi3_0/TrigTestDeltaPhiGun.root","DeltaPhiGunPt100dR3_0",10.,true));
//	wrapperVect.push_back(new AnalysisWrapper("deltaPhiGun/DeltaPhiGunPt100dPhi0_3/TrigTestDeltaPhiGun.root","DeltaPhiGunPt100dR0_3",10.,true));
//	wrapperVect.push_back(new AnalysisWrapper("deltaPhiGun/DeltaPhiGunPt100dPhi0_1/TrigTestDeltaPhiGun.root","DeltaPhiGunPt100dR0_1",10.,true));
//	wrapperVect.push_back(new AnalysisWrapper("deltaPhiGun/DeltaPhiGunPt100dPhi0_05/TrigTestDeltaPhiGun.root","DeltaPhiGunPt100dR0_05",10.,true));
//	wrapperVect.push_back(new AnalysisWrapper("deltaPhiGun/DeltaPhiGunPt100dPhi0_01/TrigTestDeltaPhiGun.root","DeltaPhiGunPt100dR0_01",10.,true));
//	wrapperVect.push_back(new AnalysisWrapper("deltaPhiGun/DeltaPhiGunPt100dPhi0_005/TrigTestDeltaPhiGun.root","DeltaPhiGunPt100dR0_005",10.,true));

	wrapperVect[0]->analyseBti();
	wrapperVect[0]->analyseTraco();
	wrapperVect[0]->analyseGenParticles();
	wrapperVect[0]->savePlots();

	//Create a vector with colors to access during scenario iterations
	std::vector<int> colorVector;
	colorVector.push_back(kBlack);
	colorVector.push_back(kBlue);
	colorVector.push_back(kViolet);
	colorVector.push_back(kRed);
	colorVector.push_back(kOrange);
	colorVector.push_back(kGreen);
	colorVector.push_back(kTeal);
	colorVector.push_back(kCyan - 1);

	//Create a vector with scenarios to access during scenario iterations
	std::vector<TString> legendNameVector;
	legendNameVector.push_back(TString("#Delta#Phi 3.0"));
	legendNameVector.push_back(TString("#Delta#Phi 0.3"));
	legendNameVector.push_back(TString("#Delta#Phi 0.1"));
	legendNameVector.push_back(TString("#Delta#Phi 0.05"));
	legendNameVector.push_back(TString("#Delta#Phi 0.01"));
	legendNameVector.push_back(TString("#Delta#Phi 0.005"));

	std::vector<std::vector<TH1*> > btiScenariosBestCase;
	std::vector<std::vector<TH1*> > btiScenariosHtrg;

	for(unsigned int i = 0 ; i < wrapperVect.size(); i++ ){
		btiScenariosBestCase.push_back(wrapperVect[i]->analyseBtiTrigPerStatAndSlBestCase());
		btiScenariosHtrg.push_back(wrapperVect[i]->analyseBtiTrigPerStatAndSlHtrg());
	}

	TCanvas* btiBestCaseCanvas = cManager->getDividedCanvas(1,1);
	btiBestCaseCanvas->SetName("btiTrigsPerStatAndSlTogetherBestCase");
	TCanvas* btiHtrgCanvas = cManager->getDividedCanvas(1,1);
	btiHtrgCanvas->SetName("btiTrigsPerStatAndSlTogetherHtrg");

	//Get the size of the first vector stored in the scenarios
	for (int i = 0 ; i < btiScenariosBestCase[0].size()/8. ; i++){
		int scenarioCounter = 0;
		TLegend* lBestCase = new TLegend(0.5,.525,0.7,0.9);
		TLegend* lBestCase2 = new TLegend(0.7,.525,0.9,0.9);
		TLegend* lHtrg = new TLegend(0.5,.525,0.7,0.9);
		TLegend* lHtrg2 = new TLegend(0.7,.525,0.9,0.9);

		btiBestCaseCanvas->cd(i+1)->SetLogy();
		//Build the frame
		TString frameName = "BTI Trigger (Best Case) f#ddot{u}r Station ";
		frameName += i+1;
		frameName += ";# BTI Trigger;Relative H#ddot{a}ufigkeit";

		TH2D* frame = new TH2D("frameBtiBestCase",frameName,1,-0.5,10.5,1,0.000001,1);
		frame->SetStats(kFALSE);
		frame->Draw();

		for (std::vector<std::vector<TH1*> >::const_iterator scenarioIt = btiScenariosBestCase.begin();
				scenarioIt != btiScenariosBestCase.end() ; scenarioIt++){

			//Loop over the scenarios
			scenarioIt->at(2*i)->SetLineColor(colorVector[scenarioCounter]);
			scenarioIt->at(2*i)->Scale(1/scenarioIt->at(2*i)->Integral());
			scenarioIt->at(2*i)->Draw("same");
			scenarioIt->at(2*i+1)->SetLineColor(colorVector[scenarioCounter]);
			scenarioIt->at(2*i+1)->SetLineStyle(2);
			scenarioIt->at(2*i+1)->Scale(1/scenarioIt->at(2*i+1)->Integral());
			scenarioIt->at(2*i+1)->Draw("same");

			lBestCase->AddEntry(scenarioIt->at(2*i),(legendNameVector[scenarioCounter] + " SL 1").Data());
			lBestCase2->AddEntry(scenarioIt->at(2*i+1),(legendNameVector[scenarioCounter] + " SL 3").Data());

			scenarioCounter++;
		}
		lBestCase->Draw();
		lBestCase2->Draw();

		btiBestCaseCanvas->Update();
		TString plotTitle("BTI Trigger (Best Case) f#ddot{u}r Station ");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    gr1_band[jj][nbin]->SetPoint(52,0,1);
    gr1_band[jj][nbin]->SetFillColor(color[nbin]-10);
  }
}


  ofstream ff("r2_eta_pPb.txt");

  for(int nbin=1;nbin<2;nbin++) {
    cout<<"i="<<nbin+1<<endl;
    ff<<"i="<<nbin+1<<endl;
    for(int j=0;j<5;j++) { cout<<slope[nbin][j]<<" "; ff<<slope[nbin][j]<<", ";}
    cout<<endl;
    ff<<endl;
    for(int j=0;j<5;j++) { cout<<slope_err[nbin][j]<<" "; ff<<slope_err[nbin][j]<<", ";}
    cout<<endl;
    ff<<endl;
  }

  TString histtitle[5] = {"120 #leq N_{trk}^{offline} < 150", "150 #leq N_{trk}^{offline} < 185", "185 #leq N_{trk}^{offline} < 220", "220 #leq N_{trk}^{offline} < 260", "N_{trk}^{offline} #geq 260"};

  TLine* ll = new TLine(0,1.,2.65,1.);
  TCanvas* c = new TCanvas("c","c",530,500);
  makeMultiPanelCanvas(c,2,2,0.01,0.0,0.2,0.18,0.02);
  TH2D* htmp = new TH2D("htmp",";#eta^{a};#sqrt{r_{2}(-#eta^{a},-#eta^{b}) #times r_{2}(#eta^{a},#eta^{b})}",100,-0.1,2.5,100,0.601,1.06);
  fixedFontHist(htmp,1.8,2.1);
  htmp->GetXaxis()->CenterTitle();
  htmp->GetYaxis()->CenterTitle();
  htmp->GetYaxis()->SetTitleSize(htmp->GetYaxis()->GetTitleSize()*1.2);
  htmp->GetXaxis()->SetTitleSize(htmp->GetXaxis()->GetTitleSize()*1.);
  for(int jj=0;jj<4;jj++)
  {
    c->cd(jj+1);
    htmp->Draw();
    ll->Draw("Lsame");

    gr_tot[jj][1]->SetMarkerStyle(25);
    gr_tot[jj][1]->Draw("Psame");
    fit[jj][1]->Draw("Lsame");

    gr_tot[jj+5][1]->SetMarkerStyle(24);
    gr_tot[jj+5][1]->Draw("Psame");

    gr_tot[jj+10][1]->SetMarkerStyle(28);
    gr_tot[jj+10][1]->SetMarkerColor(1);
    gr_tot[jj+10][1]->Draw("Psame");

/*
    gr_tot[jj+8][1]->SetMarkerStyle(28);
    fit[jj+8][1]->SetLineStyle(9);
    gr_tot[jj+8][1]->Draw("Psame");
    fit[jj+8][1]->Draw("Lsame");
*/
  }
  TLatex* latex = new TLatex();
  latex->SetNDC();
  latex->SetTextSize(1.3*latex->GetTextSize());
  c->cd(1);
  latex->DrawLatex(0.27,0.05,histtitle[0]);
//  latex->DrawLatex(0.58,0.74,"0.3 < p_{T}^{a} < 3 GeV/c");
//  latex->DrawLatex(0.58,0.65,"p_{T}^{b} > 0 GeV/c");
  latex->DrawLatex(0.24,0.885,"CMS pPb #sqrt{s_{NN}} = 5.02 TeV");
  latex->DrawLatex(0.68,0.7,"L_{int} = 35 nb^{-1}");
  c->cd(2);
  latex->DrawLatex(0.12,0.27,"0.3 < p_{T}^{a} < 3 GeV/c");
  latex->DrawLatex(0.12,0.18,"p_{T}^{b} > 0 GeV/c");

//.........这里部分代码省略.........
            }
        }
        else {
            if(LastSeedNumber != -1) {
                if(LastSimTrackvalid == 1 && LastPassSegmentFilter == true) {
                    if(debug) cout << "Filling valid track efficiency " << LastSeedNumber << endl;
                    SeedNumberforSimTrackvalidHist->Fill(LastSeedNumber);
                    SeedEfficiencyforSimTrackvalidHist->Fill(LastPurityFull==true?1:0);
                }
                else {
                    SeedNumberforSimTrackinvalidHist->Fill(LastSeedNumber);
                    SeedEfficiencyforSimTrackinvalidHist->Fill(LastPurityFull==true?1:0);
                }
                LastPurityFull = false;
            }
        }
        LastPassSegmentFilter = PassSegmentFilter;
        LastSeedNumber = SeedNumber;
        LastSimTrackvalid = SimTrackvalid;
    }
    if(LastSeedNumber != -1) {
        if(LastSimTrackvalid == 1 && LastPassSegmentFilter == true) {
            SeedNumberforSimTrackvalidHist->Fill(LastSeedNumber);
            SeedEfficiencyforSimTrackvalidHist->Fill(LastPurityFull==true?1:0);
        }
        else {
            SeedNumberforSimTrackinvalidHist->Fill(LastSeedNumber);
            SeedEfficiencyforSimTrackinvalidHist->Fill(LastPurityFull == true?1:0);
        }
    }

    TCanvas* SimTrackvalidCanvas = new TCanvas("SimTrackvalidCanvas", "SimTrackvalidCanvas", 800, 600);
    SimTrackvalidCanvas->cd();
    SimTrackvalidHist->Draw();
    string SimTrackvalidCanvasName = FinalOutput + "SimTrackvalid" + OutputPlotNameFix;
    SimTrackvalidCanvas->SaveAs(SimTrackvalidCanvasName.c_str());

    TCanvas* SeedPtforSimTrackvalidCanvas = new TCanvas("SeedPtforSimTrackvalidCanvas", "SeedPtforSimTrackvalidCanvas", 800, 600);
    SeedPtforSimTrackvalidCanvas->cd();
    SeedPtforSimTrackvalidHist->Draw();
    string SeedPtforSimTrackvalidCanvasName = FinalOutput + "SeedPtforSimTrackvalid" + OutputPlotNameFix;
    SeedPtforSimTrackvalidCanvas->SaveAs(SeedPtforSimTrackvalidCanvasName.c_str());

    TCanvas* SeeddeltaPtforSimTrackvalidCanvas = new TCanvas("SeeddeltaPtforSimTrackvalidCanvas", "SeeddeltaPtforSimTrackvalidCanvas", 800, 600);
    SeeddeltaPtforSimTrackvalidCanvas->cd();
    SeeddeltaPtforSimTrackvalidHist->SetStats(1);
    gStyle->SetOptFit(0111);
    //SeeddeltaPtforSimTrackvalidHist->Fit("gaus", "", "", -1., 1.);
    SeeddeltaPtforSimTrackvalidHist->GetXaxis()->SetTitle("(recPt-simPt)/simPt");
    SeeddeltaPtforSimTrackvalidHist->GetXaxis()->CenterTitle();
    SeeddeltaPtforSimTrackvalidHist->Draw();
    string SeeddeltaPtforSimTrackvalidCanvasName = FinalOutput + "SeeddeltaPtforSimTrackvalid" + OutputPlotNameFix;
    SeeddeltaPtforSimTrackvalidCanvas->SaveAs(SeeddeltaPtforSimTrackvalidCanvasName.c_str());

    TCanvas* SeeddeltaPhiforSimTrackvalidCanvas = new TCanvas("SeeddeltaPhiforSimTrackvalidCanvas", "SeeddeltaPhiforSimTrackvalidCanvas", 800, 600);
    SeeddeltaPhiforSimTrackvalidCanvas->cd();
    SeeddeltaPhiforSimTrackvalidHist->SetStats(1);
    SeeddeltaPhiforSimTrackvalidHist->GetXaxis()->SetTitle("(recPhi-simPhi)");
    SeeddeltaPhiforSimTrackvalidHist->GetXaxis()->CenterTitle();
    SeeddeltaPhiforSimTrackvalidHist->Draw();
    string SeeddeltaPhiforSimTrackvalidCanvasName = FinalOutput + "SeeddeltaPhiforSimTrackvalid" + OutputPlotNameFix;
    SeeddeltaPhiforSimTrackvalidCanvas->SaveAs(SeeddeltaPhiforSimTrackvalidCanvasName.c_str());

    TCanvas* SeeddeltaEtaforSimTrackvalidCanvas = new TCanvas("SeeddeltaEtaforSimTrackvalidCanvas", "SeeddeltaEtaforSimTrackvalidCanvas", 800, 600);
    SeeddeltaEtaforSimTrackvalidCanvas->cd();
    SeeddeltaEtaforSimTrackvalidHist->SetStats(1);

//.........这里部分代码省略.........
    pairDeltaEtaDeltaPhiHistogramNormalized = (TH2D*)toyV2GeneratorFile->Get("pairDeltaEtaDeltaPhiHistogramNormalized");
    if(!pairDeltaEtaDeltaPhiHistogramNormalized) {
      cerr << "\n Cannot find pairDeltaEtaDeltaPhiHistogramNormalized" << endl;
      return;
    }

  } // check on readHistogramFile

  TCanvas *c1 = new TCanvas("c1", "Plotting Pair Events", 200, 10, 700, 500);
  //
  // Use ROOT's gStyle function to set up plot options
  //
  gStyle->SetOptStat(1110);
  gStyle->SetStatTextColor(kBlue);
  gStyle->SetFuncColor(kBlue);
  gStyle->SetFuncWidth(2);
  gStyle->SetOptFit(1111);
  gStyle->SetFitFormat("6.5g");
	
  c1->SetFillColor(kWhite);

  c1->Divide(2,2);

  c1->cd(1);
  
  if(useFlow) {
    particlePhiHistogram->SetStats(0);
    particlePhiHistogram->SetMinimum(0);
    float fMaximum = particlePhiHistogram->GetMaximum();
    int iMaximum = 1.5*fMaximum;
    particlePhiHistogram->SetMaximum(iMaximum);
		
    particlePhiHistogram->SetLineColor(4);
    particlePhiHistogram->Draw();
    particlePhiLabHistogram->SetLineColor(1);
    particlePhiLabHistogram->Draw("same");
		
    char labelBufferTitle[200];
    sprintf(labelBufferTitle, "Events = %d, Particles/Event = %d, V2 factor = %4.2f", numberEvents, particlesPerEvent, v2Factor);
    TLatex *text1 = new TLatex(0.5, 0.90*iMaximum, labelBufferTitle);
    text1->SetTextColor(4);
    text1->SetTextSize(0.045);
    text1->Draw();
		
    char labelBuffer2[200];
    sprintf(labelBuffer2, "Blue: azimuthal distribution relative to reaction plane");
    TLatex *text2 = new TLatex(0.5, 0.20*iMaximum, labelBuffer2);
    text2->SetTextColor(4);
    text2->SetTextSize(0.035);
    text2->Draw();
		
    char labelBuffer3[200];
    sprintf(labelBuffer3, "Black: azimuthal distribution relative to lab frame");
    TLatex *text3 = new TLatex(0.5, 0.10*iMaximum, labelBuffer3);
    text3->SetTextColor(1);
    text3->SetTextSize(0.035);
    text3->Draw();
		
  }
  else {
    particlePhiLabHistogram->SetStats(0);
    float fMaximum = particlePhiLabHistogram->GetMaximum();
    int iMaximum = 1.5*fMaximum;
    particlePhiLabHistogram->SetMaximum(iMaximum);
    particlePhiLabHistogram->SetMinimum(0);
    particlePhiLabHistogram->Draw();

void clusterGaussReso() {

  TFile* file = new TFile("complete.root");

  TTree* tree = (TTree *)file->Get("Signal");

  int nEvents = tree->GetEntries();

  double addresses[121] = {};
  for (int k=0; k<121; k++){
    std::stringstream ss2;
    ss2 << k; 
    string str = "Crystal_"+ss2.str();
    const char* charstr = str.c_str(); 
    tree->SetBranchAddress(charstr, &addresses[k]);
  }

  TH1D* energyResoG = new TH1D("energyResoG", "Gauss_Energy_Resolution",
			      100, -10, 10);

  TH2D* posResoG = new TH2D("posResoG", "Gauss_Pos_Resolution", 
			    11, -5, 5, 11, -5, 5);

  double energyThresHi = 50.; //set energy threshold to start looking for bumps
  double energyThresLo = .5; //energy lower bound


  //Going through each scan of the calorimeter
  for (int k=0; k<nEvents; k++)
    {
      tree->GetEntry(k);

      pair<int, double> bump(0, 0.); 
      vector<pair<int, double> > geant;
      vector<pair<int, double> > hitMap;

      for(int i=0; i<121; i++)
	{
	  pair<int, double> hit(i, addresses[i]);
	  geant.push_back(hit);
	  if (addresses[i] > energyThresHi)
	    { hitMap.push_back(hit);}
	  if (hit.second>bump.second)
	    { bump = hit;}
	}

      gaussFit(bump, geant, energyResoG, posResoG);


    } //end of event

 TCanvas* canvas = new TCanvas("canvas", "canvas", 700,700);
 canvas->Divide(2, 1);
 TF1* g1d = new TF1("g1d", "gaus", -10, 10);
 TF2* g2d2 = new TF2("g2d2", Gaus2D, -10, 10, -10, 10, 6);
 g2d2->SetParameters(0., 0., 1., 1., 0., 100000);
 g2d2->FixParameter(4, 0.);
 energyResoG->Fit(g1d);
canvas->cd(1); energyResoG->Draw();
 posResoG->Fit(g2d2);
canvas->cd(2); posResoG->Draw("BOX");

}

void plot_eff_NoErrors(const string& fInputFile, const string& fPlot, const double fOP, const string& fTitle, const string& fXAxisTitle, const string& fYAxisTitle, const string& fOutputFile,
                       const double fYmin=0., const double fYmax=1.0, const Int_t fLogy=0, const Double_t fTitleOffsetX=1.0, const Double_t fTitleOffsetY=1.0,
                       const Double_t fLeftMargin=0.12, const Double_t fTopMargin=0.07, const Double_t fPlotWidth=0.8)
{
  gROOT->SetBatch(kTRUE);
  setEXOStyle();
  gStyle->SetGridColor(kGray);
  gStyle->SetOptStat(kFALSE);
  gStyle->SetPadTopMargin(fTopMargin);
  gStyle->SetPadBottomMargin(1.-fTopMargin-fPlotWidth);
  gStyle->SetPadLeftMargin(fLeftMargin);
  gStyle->SetPadRightMargin(1.-fLeftMargin-fPlotWidth);
  gROOT->ForceStyle();

  TFile *file = new TFile(fInputFile.c_str());

  TH1D *h1_total = new TH1D("h1_total","h1_total",4,0.,0.8);
  TH1D *h1_subjet1 = new TH1D("h1_subjet1","h1_subjet1",4,0.,0.8);
  TH1D *h1_subjet2 = new TH1D("h1_subjet2","h1_subjet2",4,0.,0.8);
  TH1D *h1_subjet12 = new TH1D("h1_subjet12","h1_subjet12",4,0.,0.8);

  string bin_strings[4] = {"0to0p2", "0p2to0p4", "0p4to0p6", "0p6to0p8"};

  for(int i=0; i<4; ++i)
  {
     TH2D *h2 = (TH2D*)file->Get((fPlot + bin_strings[i]).c_str());

     int bin = h2->GetXaxis()->FindBin(fOP);

     h1_total->SetBinContent(   i+1,h2->Integral(0,101,0,101));
     h1_subjet1->SetBinContent( i+1,h2->Integral(bin,101,0,101));
     h1_subjet2->SetBinContent( i+1,h2->Integral(0,101,bin,101));
     h1_subjet12->SetBinContent(i+1,h2->Integral(bin,101,bin,101));
  }

  TCanvas *c = new TCanvas("c", "",1000,800);
  c->cd();

  TH2D *bkg = new TH2D("bkg","",100,0.,0.8,100,fYmin,fYmax);
  bkg->GetXaxis()->SetTitle(fXAxisTitle.c_str());
  bkg->GetYaxis()->SetTitle(fYAxisTitle.c_str());
  bkg->SetTitleOffset(fTitleOffsetX,"X");
  bkg->SetTitleOffset(fTitleOffsetY,"Y");
  bkg->Draw();
  c->SetGridx();
  c->SetGridy();

  TGraph *g_eff_subjet1 = new TGraph(4);
  g_eff_subjet1->SetMarkerStyle(20);
  g_eff_subjet1->SetMarkerColor(kBlue+1);
  g_eff_subjet1->SetLineColor(kBlue+1);
  g_eff_subjet1->SetLineStyle(2);
  for(int i=0; i<4; ++i)
    g_eff_subjet1->SetPoint(i,h1_total->GetBinCenter(i+1),h1_subjet1->GetBinContent(i+1)/h1_total->GetBinContent(i+1));

  TGraph *g_eff_subjet2 = new TGraph(4);
  g_eff_subjet2->SetMarkerStyle(21);
  g_eff_subjet2->SetMarkerColor(kGreen+1);
  g_eff_subjet2->SetLineColor(kGreen+1);
  g_eff_subjet2->SetLineStyle(2);
  for(int i=0; i<4; ++i)
    g_eff_subjet2->SetPoint(i,h1_total->GetBinCenter(i+1),h1_subjet2->GetBinContent(i+1)/h1_total->GetBinContent(i+1));

  TGraph *g_eff_subjet12 = new TGraphAsymmErrors(4);
  g_eff_subjet12->SetMarkerStyle(22);
  g_eff_subjet12->SetMarkerColor(kBlack);
  g_eff_subjet12->SetLineColor(kBlack);
  g_eff_subjet12->SetLineStyle(2);
  for(int i=0; i<4; ++i)
    g_eff_subjet12->SetPoint(i,h1_total->GetBinCenter(i+1),h1_subjet12->GetBinContent(i+1)/h1_total->GetBinContent(i+1));

  TGraph *g_eff_subjet12_prod = new TGraph(4);
  g_eff_subjet12_prod->SetMarkerStyle(26);
  g_eff_subjet12_prod->SetMarkerColor(kRed);
  g_eff_subjet12_prod->SetLineColor(kRed);
  g_eff_subjet12_prod->SetLineStyle(2);
  for(int i=0; i<4; ++i)
    g_eff_subjet12_prod->SetPoint(i,h1_total->GetBinCenter(i+1),(h1_subjet1->GetBinContent(i+1)*h1_subjet2->GetBinContent(i+1))/(h1_total->GetBinContent(i+1)*h1_total->GetBinContent(i+1)));

  g_eff_subjet1->Draw("PLsame");
  g_eff_subjet2->Draw("PLsame");
  g_eff_subjet12->Draw("PLsame");
  g_eff_subjet12_prod->Draw("PLsame");

  TLegend *legend = new TLegend(.6,.3,.85,.55);
  legend->SetBorderSize(0);
  legend->SetFillColor(0);
  legend->SetFillStyle(0);
  legend->SetTextFont(42);
  legend->SetTextSize(0.05);
  legend->AddEntry(g_eff_subjet1, "#varepsilon_{1}","lp");
  legend->AddEntry(g_eff_subjet2, "#varepsilon_{2}","lp");
  legend->AddEntry(g_eff_subjet12, "#varepsilon_{1 & 2}","lp");
  legend->AddEntry(g_eff_subjet12_prod, "#varepsilon_{1} #times #varepsilon_{2}","lp");
  legend->Draw();

  TLatex l1;
  l1.SetTextAlign(13);
  l1.SetTextFont(42);
  l1.SetTextSize(0.045);
//.........这里部分代码省略.........

float fitAndDraw( const std::string& outputdir, TH1F* h1 ) {


  float returnConst=0.;

  TCanvas* c1 = new TCanvas("c1", "", 600, 600);
  c1->cd();

  float xMin = h1->GetXaxis()->GetXmin();
  float xMax = h1->GetXaxis()->GetXmax();
  float yMax = 1.1*h1->GetMaximum();


  TH2D* axes = new TH2D( "axes", "", 10, xMin, xMax, 10, 0., yMax );
  axes->SetXTitle("Cluster Position [mm]");
  //  axes->SetXTitle(h1->GetName());
  axes->SetYTitle("Events");
  axes->Draw();

  h1->SetFillColor(29);
  h1->Draw("same");

  int opt=2;

  if( opt==1 ) { // gaussian fit

    float maxPos = h1->GetBinCenter(h1->GetMaximumBin());

    TF1* f1 = new TF1(Form("f1_%s", h1->GetName()), "gaus", -10., 10.);
    
    f1->SetParameter( 1, maxPos );
    h1->Fit(f1, "QRN");
    
    for( int i=0; i<4; ++i ) {

      float m = f1->GetParameter(1);
      float s = f1->GetParameter(2);
      float nSigmas = 1.2;
      f1->SetRange( m-nSigmas*s, m+nSigmas*s );
      if( i==3 ) 
        h1->Fit(f1, "RQ");
      else
        h1->Fit(f1, "RNQ");

    }


    f1->SetLineColor(kRed);
    returnConst = f1->GetParameter(1);

  } else if( opt==2 ) { // average of high bins

    float maximum = h1->GetMaximum();
    float thresh = 0.25*maximum;

    float total=0.;
    float denom = 0.;
    
    int nBins = h1->GetNbinsX();
    TH1F* h1_usedBins = new TH1F( Form("usedBins_%s", h1->GetName()), "", nBins, xMin, xMax );

    for( int ibin=1; ibin<nBins; ++ibin ) {

      if( h1->GetBinContent(ibin)<thresh ) continue;

      h1_usedBins->SetBinContent( ibin, h1->GetBinContent(ibin) );
  
      // average:
      total += h1->GetBinCenter(ibin)*h1->GetBinContent(ibin);
      denom += h1->GetBinContent(ibin);

      //// weighted average:
      //total += h1->GetBinCenter(ibin)*h1->GetBinContent(ibin);
      //denom += h1->GetBinContent(ibin);

    }

    TLine* lineThresh = new TLine( xMin, thresh, xMax, thresh );
    lineThresh->SetLineStyle(2);
    lineThresh->Draw("same");

    returnConst = total/denom;
    h1_usedBins->SetFillColor(kOrange);
    h1_usedBins->Draw("same");


  } else if( opt==3 ) { // maximum

    float maxBinCenter = h1->GetBinCenter(h1->GetMaximumBin());

    returnConst = maxBinCenter;

  }


  TLine* lineOffset = new TLine( returnConst, 0., returnConst, yMax );
  lineOffset->SetLineColor(kRed);
  lineOffset->SetLineWidth(3);
  lineOffset->Draw("same");

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

void run() {


  TGraphErrors* gr0 = new TGraphErrors(0);
  TGraphErrors* gr1 = new TGraphErrors(0);
  TGraphErrors* gr2 = new TGraphErrors(0);
  TGraphErrors* gr3 = new TGraphErrors(0);

  TGraphErrors* gr23 = new TGraphErrors(0);

  TGraphErrors* gr_cef3 = new TGraphErrors(0);
  gr_cef3->SetName("cef3");

  std::vector< std::pair<int, float> > run_pos;
  run_pos.push_back( std::pair< int, float>(225, 250.0) );
  run_pos.push_back( std::pair< int, float>(226, 251.0) );
  run_pos.push_back( std::pair< int, float>(227, 252.0) );
  run_pos.push_back( std::pair< int, float>(228, 253.0) );
  run_pos.push_back( std::pair< int, float>(229, 254.0) );
  run_pos.push_back( std::pair< int, float>(230, 248.0) );
  run_pos.push_back( std::pair< int, float>(231, 247.0) );
  run_pos.push_back( std::pair< int, float>(232, 246.0) );
  run_pos.push_back( std::pair< int, float>(233, 245.0) );
  run_pos.push_back( std::pair< int, float>(234, 244.0) );
  run_pos.push_back( std::pair< int, float>(235, 243.0) );
  run_pos.push_back( std::pair< int, float>(236, 242.0) );
  run_pos.push_back( std::pair< int, float>(237, 249.0) );
  run_pos.push_back( std::pair< int, float>(238, 255.0) );
  run_pos.push_back( std::pair< int, float>(239, 256.0) );
  run_pos.push_back( std::pair< int, float>(240, 257.0) );
  run_pos.push_back( std::pair< int, float>(241, 258.0) );
  run_pos.push_back( std::pair< int, float>(242, 259.0) );
  run_pos.push_back( std::pair< int, float>(243, 260.0) );
  run_pos.push_back( std::pair< int, float>(246, 261.0) );
  run_pos.push_back( std::pair< int, float>(247, 262.0) );
  run_pos.push_back( std::pair< int, float>(248, 263.0) );
  run_pos.push_back( std::pair< int, float>(249, 264.0) );
  run_pos.push_back( std::pair< int, float>(250, 265.0) );
  run_pos.push_back( std::pair< int, float>(251, 266.0) );
  run_pos.push_back( std::pair< int, float>(252, 267.0) );
  run_pos.push_back( std::pair< int, float>(253, 268.0) );
  run_pos.push_back( std::pair< int, float>(254, 270.0) );
  //run_pos.push_back( std::pair< int, float>(255, 272.5) );
  run_pos.push_back( std::pair< int, float>(244, 275.0) );

  for( unsigned i=0; i<run_pos.size(); ++i ) {

    hodoScan p(run_pos[i].first);
    std::vector<float> eff = p.Loop();
    gr0->SetPoint( i, run_pos[i].second, eff[0] );
    gr1->SetPoint( i, run_pos[i].second, eff[1] );
    gr2->SetPoint( i, run_pos[i].second, eff[2] );
    gr3->SetPoint( i, run_pos[i].second, eff[3] );
    
    gr23->SetPoint( i, run_pos[i].second, 0.5*(eff[3]+eff[2]) );
    
    cef3Scan c(run_pos[i].first);
    TH1F* h1_cef3 = c.Loop();
    
    //TF1* f1 = new TF1(Form("f1_%d", run_pos[i].first), "gaus", 3000., 6000.);
    //h1_cef3->Fit( f1, "RQL" );
    //for( unsigned j=0; j<4; ++j ) {
    //  float mu = f1->GetParameter(1);
    //  float sigma = f1->GetParameter(2);
    //  f1->SetRange( mu-sigma, mu+sigma );
    //  h1_cef3->Fit( f1, "RQL" );
    //}

    //TCanvas* c2 = new TCanvas("c2", "", 600, 600);
    //c2->cd();
    //h1_cef3->Draw();
    //c2->SaveAs(Form("provaY/%d.png", run_pos[i].first));
    //delete c2;
    //gr_cef3->SetPoint( i, run_pos[i].second, f1->GetParameter(1) );
    //gr_cef3->SetPointError( i, 0., f1->GetParError(1) );
  

    gr_cef3->SetPoint( i, run_pos[i].second, h1_cef3->GetBinCenter(h1_cef3->GetMaximumBin()));
    
  }

  TCanvas* c1 = new TCanvas("c1", "", 600, 600);
  c1->cd();

  TH2D* axes = new TH2D("axes", "", 10, 200., 300., 10, 0., 0.15 );
  axes->SetXTitle("Y Position [mm]");
  axes->SetYTitle("Efficiency");
  axes->Draw();

  gr0->SetMarkerStyle(20);
  gr0->SetMarkerColor(46);
  gr0->SetMarkerSize(1.);

  gr1->SetMarkerStyle(21);
  gr1->SetMarkerColor(29);
  gr1->SetMarkerSize(1.);

  gr2->SetMarkerStyle(22);
  gr2->SetMarkerColor(38);
  gr2->SetMarkerSize(1.);
//.........这里部分代码省略.........