C++ TGraph::GetY方法代码示例

void plot_mc(string name="g", double q2=1.9, Int_t n=0)
{
//    SetAtlasStyle();
    TGraphErrors *p;

    TString nn;
    nn.Form("ave_%s_vs_x_for_Q2=%g",name.c_str(),q2);
    gDirectory->GetObject(nn,p);

//    p->Print();

    Double_t ratsize = 0.0;

    TCanvas *c = new TCanvas("PDF","pdf",600,600);
    TPad *pad1 = new TPad("pad1","pad1",0.,ratsize,1.,1.);
    pad1->SetLogx();
    pad1->Draw();
    pad1->cd();
    p->GetXaxis()->Set(101,0.0001,1.);
    p->SetFillColor(kRed-2);
    p->SetFillStyle(3001);
    p->SetLineWidth(1);
    p->SetLineColor(kRed);
    p->GetYaxis()->SetTitle(name.c_str());

    p->GetYaxis()->SetTitleSize(0.06);
    p->GetYaxis()->SetTitleOffset(1.);
    p->GetXaxis()->Set(101,0.0001,1.);
    p->Draw("ALE3");

    Double_t av = 0;
    Double_t av2 = 0;

    for (Int_t i = 1; i<n+1 ; i++) {
        nn.Form("%s_vs_x_for_Q^{2}=%g_%i",name.c_str(),q2,i);
        TGraph *pp = NULL;
        gDirectory->GetObject(nn,pp);
        if (pp != NULL) {
            pp->SetLineColor(kGreen);
            pp->Draw("L");

            av  += pp->GetY()[0];
            av2 += pp->GetY()[0]*pp->GetY()[0];
            cout << i << " "<<pp->GetY()[0] << endl;


        }
    }

    av  /= n;
    av2 = sqrt(av2/n - av*av);
//    cout << n << " "<<av << " "<< av2<<endl;


    p->Draw("E3C");


}

void
plotTanb(const char* filename, const char* channel, bool draw_injected_=false, double min_=0., double max_=60., bool log_=false, std::string dataset_="CMS Preliminary,  H#rightarrow#tau#tau,  4.9 fb^{-1} at 7 TeV, 19.8 fb^{-1} at 8 TeV", std::string xaxis_="m_{A} [GeV]", std::string yaxis_="#bf{tan#beta}")
{
  TFile* file = TFile::Open(filename);
  // retrieve TGraphs from file
  TGraph* expected = (TGraph*)file->Get(std::string(channel).append("/expected").c_str());
  TGraph* expected_low = (TGraph*)file->Get(std::string(channel).append("/expected_low").c_str());
  TGraph* observed = (TGraph*)file->Get(std::string(channel).append("/observed").c_str());
  TGraph* observed_low = (TGraph*)file->Get(std::string(channel).append("/observed_low").c_str());
  TGraph* upperLEP = (TGraph*)file->Get(std::string(channel).append("/upperLEP").c_str());
  TGraph* lowerLEP = (TGraph*)file->Get(std::string(channel).append("/lowerLEP").c_str());
  TGraphAsymmErrors* innerBand = (TGraphAsymmErrors*)file->Get(std::string(channel).append("/innerBand").c_str());
  TGraphAsymmErrors* innerBand_low = (TGraphAsymmErrors*)file->Get(std::string(channel).append("/innerBand_low").c_str());
  TGraphAsymmErrors* outerBand = (TGraphAsymmErrors*)file->Get(std::string(channel).append("/outerBand").c_str());
  TGraphAsymmErrors* outerBand_low = (TGraphAsymmErrors*)file->Get(std::string(channel).append("/outerBand_low").c_str()); 
  TGraphAsymmErrors* plain = (TGraphAsymmErrors*)file->Get(std::string(channel).append("/plain").c_str());
  TGraphAsymmErrors* plain_low = (TGraphAsymmErrors*)file->Get(std::string(channel).append("/plain_low").c_str());

  // this is new for injected plot 
  TGraph* injected = 0;;
  if(draw_injected_) {injected = (TGraph*)file->Get("injected/observed");}

  if(draw_injected_){
    for( int i=0; i<expected->GetN(); i++){
      double shift = injected->GetY()[i]-expected->GetY()[i];
      innerBand->SetPoint(i, innerBand->GetX()[i], innerBand->GetY()[i]+shift);
      outerBand->SetPoint(i, outerBand->GetX()[i], outerBand->GetY()[i]+shift);
    }
  }

  // this functionality is not yet supported
  std::map<double, TGraphAsymmErrors*> higgsBands;
  // this functionality is not yet supported
  std::map<std::string, TGraph*> comparisons;

  // set up styles
  SetStyle();
  // do the plotting 
  TCanvas canv = TCanvas("canv", "Limits", 600, 600);
  // do the plotting 
  plottingTanb(canv, plain, plain_low, innerBand, innerBand_low, outerBand, outerBand_low, expected, expected_low, observed, observed_low, lowerLEP, upperLEP, higgsBands, comparisons, xaxis_, yaxis_, injected, min_, max_, log_);
  /// setup the CMS Preliminary
  CMSPrelim(dataset_.c_str(), "", 0.145, 0.835);
  // write results to files
  canv.Print(std::string(channel).append("_tanb").append(".png").c_str());
  canv.Print(std::string(channel).append("_tanb").append(".pdf").c_str()); 
  canv.Print(std::string(channel).append("_tanb").append(".eps").c_str()); 
  return;
}

void BoundaryFinder::StoreDataFromGraph(const TGraph &gr)
{
    Double_t *x = gr.GetX();
    Double_t *y = gr.GetY();

    ComputeCenter(gr);

    Double_t xr=0, yr=0;

    for (UInt_t i=0; i<gr.GetN(); i++)
    {
        xr = x[i]-fCenter.first;
        yr = y[i]-fCenter.second;
        point_t p = {::sqrt(xr*xr + yr*yr), i};
        fSortedVals.insert(std::make_pair(::atan2(yr,xr),p));
    }
}

void BoundaryFinder::ComputeCenter(const TGraph& gr)
{
    Double_t *x = gr.GetX();
    Double_t *y = gr.GetY();
    UInt_t n = gr.GetN();

    Double_t sum_x = 0;
    Double_t sum_y = 0;
    for (UInt_t i=0; i<n; i++)
    {
        sum_x += x[i];
        sum_y += y[i];
    }

    fCenter.first = sum_x/n;
    fCenter.second = sum_y/n;
}

//______________________________________________________________________________
void Draw(TFile* f, const char* gname, TLegend* l, Bool_t normalized)
{
  if (!f) return;

  TGraph* g = static_cast<TGraph*>(f->Get(gname));

  if (!g) return;

  if ( normalized )
  {
    g = static_cast<TGraph*>(g->Clone());
    for ( Int_t i = 0; i < g->GetN(); ++i )
    {
      Double_t y = g->GetY()[i];
      g->SetPoint(i,g->GetX()[i],y/NTOTALNUMBEROFPADS);
    }
  }

  g->Draw("lp");
  g->GetXaxis()->SetNdivisions(505);
  g->GetXaxis()->SetNoExponent();

  if (l) l->AddEntry(g,gname,"LP");
}

//.........这里部分代码省略.........
     prettyChain.GetEntry(ientry);

     //cout << "realTime is " << header->realTime << endl; 
     //cout << " here " << endl; 
   
     //This step produces a calibrated UsefulAnitaEvent
     UsefulAnitaEvent realEvent(event,WaveCalType::kDefault,hk);
     //cout << " after useful " << endl; 

     //cout << realEvent.eventNumber << " " << header->eventNumber << endl;
     //  cout << realEvent.gotCalibTemp << " " << realEvent.calibTemp << endl;
     count++;

     //initialize 
     double min_pkpk = 0.0;
    
     //for (int j = 0; j < num_ant; j++)
     //{
     //pkpk[j] = 0.0; 
     //} 

     for (int j = 0; j < num_phi; j++)
       {
	 pkpk[j] = 0.0; 
       } 

    //for (int iant = 0; iant < num_ant; iant++)
    //{ 
    
    //TGraph *gr = new TGraph(0); 
    // gr = realEvent.getGraph(iant,pol);
    //pkpk[iant] = (gr->GetY()[TMath::LocMax(gr->GetN(),gr->GetY())]) - (gr->GetY()[TMath::LocMin(gr->GetN(),gr->GetY())]);

    //cout << pkpk[iant] << endl; 

    //delete gr;

      //} //loop over antennas
     
     max_index = 0;
     min_index = 0;         
     max = 0.0;
     min = 0.0;

     for (int iphi = 0; iphi < num_phi; iphi++)
      {
	max_index = 0;
	min_index = 0;         
	max = 0.0;
	min = 0.0;

	TGraph *gr = new TGraph(0);
	gr = realEvent.getGraph(ring,iphi,pol); 

	max_index = TMath::LocMax(gr->GetN(),gr->GetY());
	min_index = TMath::LocMin(gr->GetN(),gr->GetY());         
	max = gr->GetY()[max_index];
	min = gr->GetY()[min_index];  
	
	pkpk[iphi] = max - min;  

	//cout << pkpk[iphi] << endl; 

	//cout << max << "   " << min << endl; 
	//cout << max_index << "   " << min_index << endl; 

	delete gr; 
      } // loop over phi sectors 

    //min_pkpk = pkpk[TMath::LocMin(num_ant,pkpk)];  
    min_pkpk = pkpk[TMath::LocMin(num_phi,pkpk)]; 
    //cout << "min pkpk is " << min_pkpk << endl; 

    hmin_signal->Fill(min_pkpk);

    if (min_pkpk < 50) { cout << "min pkpk is " << min_pkpk << "run is " << header->run << "event number is " << header->eventNumber << endl; } 

  } //loop over events ends

  cerr << endl;
  cout << "Processed " << count << " events.\n";

   //just to check how a graph looks
  //TCanvas *c = new TCanvas("c","c",1000,800); 
  //gr->Draw("alp"); 
  //gr->GetXaxis()->SetTitle("Time (ns)");
  //gr->GetYaxis()->SetTitle("Voltage (mV)"); 
  //gr->SetTitle("Voltage-time waveform using UsefulAnitaEvent"); 
  //c->SaveAs("gr.png"); 
  //delete c; 

  TCanvas *h = new TCanvas("h","h",1000,800);
  h->SetLogy();  
  hmin_signal->SetStats(0); 
  hmin_signal->Draw("");
  hmin_signal->SaveAs(Form("hmin_signal_pol%iring%i.root",pol,ring)); 
  h->SaveAs(Form("min_signal_pol%iring%i.png",pol,ring));
  h->SaveAs(Form("min_signal_pol%iring%i.root",pol,ring)); 

}//end of macro

// This must run with AnitaAnalysisFramework loaded
void ALFASpectrum(int eventNumber)
{

  int run = AnitaDataset::getRunContainingEventNumber(eventNumber); 

  AnitaDataset * d = new AnitaDataset(run); 
  d->getEvent(eventNumber); 

  d->useful()->setAlfaFilterFlag(false); 

  TGraph * alfaGraph = d->useful()->getGraph(AnitaRing::kTopRing, 4, AnitaPol::kHorizontal); 

//  printf("%p\n", alfaGraph); 
  AnalysisWaveform *wf_akima = new AnalysisWaveform( alfaGraph->GetN(), alfaGraph->GetX(), alfaGraph->GetY(), 1./2.6, AnalysisWaveform::AKIMA);

//  AnalysisWaveform *wf_sparse = new AnalysisWaveform( alfaGraph->GetN(), alfaGraph->GetX(), alfaGraph->GetY(), 1./2.6, AnalysisWaveform::REGULARIZED_SPARSE_YEN);

  const TGraphAligned * uneven = wf_akima->uneven(); 

  TString fname;
  fname.Form("alfa_%d_uneven.csv",eventNumber); 

  FILE * ff = fopen(fname.Data(),"w"); 
  for (int i = 0; i < uneven->GetN(); i++)
    fprintf(ff,"%g,%g\n", uneven->GetX()[i], uneven->GetY()[i]); 
  fclose(ff); 

  fname.Form("alfa_%d_akima.csv",eventNumber); 

  const TGraphAligned * even = wf_akima->even(); 

  ff = fopen(fname.Data(),"w"); 
  for (int i = 0; i < even->GetN(); i++)
    fprintf(ff,"%g,%g\n", even->GetX()[i], even->GetY()[i]); 
  fclose(ff); 

  /*
  new TCanvas; 
  alfaGraph->Draw("al"); 
  wf_akima->setTitle("akima"); 
  wf_sparse->setTitle("sparse"); 
  wf_akima->setColor(2); 
  wf_sparse->setColor(3); 
//  wf_akima->padFreq(3); 
//  wf_sparse->padFreq(3); 
  wf_akima->drawEven("lsame"); 
  wf_sparse->drawEven("lsame"); 

  */

  TCanvas * c = new TCanvas;
  wf_akima->padEven(10); 
 //wf_sparse->padEven(9); 
//  wf_sparse->setFreqDisplayRange(0.7,0.9);
  wf_akima->drawPower("al"); 
  wf_akima->setFreqDisplayRange(0.7,0.9);
  wf_akima->setTitle(TString::Format("%d",eventNumber)); 

  TGraph *g = new TGraph(2); 
  g->GetX()[0] = 0.792; 
  g->GetX()[1] = 0.792; 
  g->GetY()[1] = 1; 
  g->GetY()[0] = 0; 


//  wf_sparse->drawPower("lsame"); 

  g->Draw("lsame"); 

  c->SaveAs(TString::Format("alfa_%d.png", eventNumber)); 


}

void PlotRakeEvolutions(const TString &sim, const TString &options="png") { 
  
#ifdef __CINT__  
  gSystem->Load("libplasma.so");
#endif

  PlasmaGlob::Initialize();
  
  // Palettes!
  gROOT->Macro("PlasmaPalettes.C");

  TString opt = options;
  
  // More makeup            
  Float_t margins[4] = {0.15,0.15,0.20,0.10};
  gStyle->SetPadLeftMargin(margins[0]);  // Margin left axis  
  gStyle->SetPadRightMargin(margins[2]);
  gStyle->SetPadTopMargin(margins[3]);  // Margin left axis  
  gStyle->SetPadBottomMargin(margins[1]);

  gStyle->SetPadTickX(1);
  gStyle->SetPadTickY(1);


  if(opt.Contains("grid")) {
    gStyle->SetPadGridX(1);
    gStyle->SetPadGridY(1);
  }

  const Int_t Nspaces = 4;
  TString phaname[Nspaces] = {"p1x1","p2x2","p3x3","x2x1"};
  TGraph *gXmean[Nspaces]; 
  TGraph *gYmean[Nspaces]; 
  TGraph *gXrms[Nspaces]; 
  TGraph *gYrms[Nspaces]; 
  TGraph *gEmit[Nspaces];
  TGraph *gCharge = NULL;
  // Special graph with an Energy spread band:
  TGraphErrors *gEneRms = NULL;
 
  Float_t maxEmit[Nspaces] = { -999., -999., -999., -999.};
  Float_t minEmit[Nspaces] = { 999., 999.,999., 999.};
  Float_t maxXmean[Nspaces] = { -999., -999., -999., -999.};
  Float_t minXmean[Nspaces] = { 999., 999., 999., 999.};
  Float_t maxXrms[Nspaces] = { -999., -999., -999., -999.};
  Float_t minXrms[Nspaces] = { 999., 999., 999., 999.};
  Float_t maxYmean[Nspaces] = { -999., -999., -999., -999.};
  Float_t minYmean[Nspaces] = { 999., 999., 999.};
  Float_t maxYrms[Nspaces] = { -999., -999., -999., -999.};
  Float_t minYrms[Nspaces] = { 999., 999., 999., 999.};
  Float_t maxCharge =-999.;
  Float_t minCharge = 999.;
  
  // Resolution:
  Int_t sizex = 600;
  Int_t sizey = 800;
  if(opt.Contains("hres")) {
    Int_t sizex = 1024;
    Int_t sizey = 768;    
  }
  

  for(Int_t i=0;i<Nspaces;i++) {
    TString filename;
    filename = Form("./%s/Plots/EmittanceEvolution/Evolutions-%s-%s.root",sim.Data(),sim.Data(),phaname[i].Data());
    
    TFile  *ifile = (TFile*) gROOT->GetListOfFiles()->FindObject(filename.Data());
    if (!ifile) ifile = new TFile(filename,"READ");
    
    
    gEmit[i] = (TGraph*) ifile->Get("gEmitvsTime");  
    gXmean[i] = (TGraph*) ifile->Get("gXmeanvsTime");
    gXrms[i] = (TGraph*) ifile->Get("gXrmsvsTime");
    gYmean[i] = (TGraph*) ifile->Get("gYmeanvsTime");
    gYrms[i] = (TGraph*) ifile->Get("gYrmsvsTime");

    // Energy spread
    if(i==0) {
      gCharge = (TGraph*) ifile->Get("gChargevsTime");  
      Int_t Npoints = gCharge->GetN();
      Double_t *yCharge = gCharge->GetY();
      for(Int_t j=0;j<Npoints;j++) {
	if(yCharge[j]>maxCharge)
	  maxCharge = yCharge[j];
	if(yCharge[j]<minCharge)
	  minCharge = yCharge[j];
      }
      Double_t *xValues = gYmean[i]->GetX();
      Double_t *yMean = gYmean[i]->GetY();
      Double_t *yRms = gYrms[i]->GetY();
      Npoints =  gYmean[i]->GetN();

      gEneRms = new TGraphErrors(Npoints,xValues,yMean,0,yRms);
      // cout << "NPOints = " << gEneRms->GetN() << endl;
      // for(Int_t j=0;j<gEneRms->GetN();j++) {
      // 	gEneRms->SetPointError(j,0.0,yRms[j]);
      // 	cout << "eoooo" << endl;
      // }
    }
    // Calculate the max and min of every set of graphs:
//.........这里部分代码省略.........

void laserCalibration(
		char* filename = "frascatirun", //input file
		int filenum = 1081, 					//file number
		int channel = 3, 						//trace channel
		int flagChannel = 5, 					//laser flag channel
		Double_t entriesN = 10,					//number of entries for prcessing
		int sleep = 10,							//sleep time between 2 processed entries, helpful for viewing traces
		bool gui = true							//enable or disable trace visualization
		)
{
	caen_5742 caen;
	Int_t nbins = 1024;
	Double_t entries = entriesN;
	Int_t bin;

	TCanvas *c1 = new TCanvas("c1","frascatirun",900,700);
	c1->Divide(1,2);
	c1->cd(1);

	TGraph* g = new TGraph();
	TH1F* lmPeaks = new TH1F("lm","Peaks Ratio", 1000, 0, 5000);
    TH1F* d = new TH1F("d","",nbins,0,nbins);
    TH1F* back = new TH1F("Back","",nbins,0,nbins);

    // input file
    char fname[100]=0;
    sprintf(fname,"%s_0%i.root",filename,filenum);

    TFile* infile = new TFile(fname);
	TTree *t = (TTree*) infile->Get("t");
	t->SetBranchAddress("caen_5742", &caen.system_clock);
	t->Print();

	if(entriesN<=0)
		entries = t->GetEntries();

	//out file
    char foutname[100]=0;
    int lm=0;
    if (channel ==3)lm=1;
    if (channel ==4)lm=2;
    sprintf(foutname,"./calibration/LM%i_out_0%i.root",lm,filenum);
	outfile = new TFile(foutname,"RECREATE");
	outTree = new TTree("LM","frascatirun output");
	calibTree = new TTree("LM_cal","frascatirun output");
	outTree->Branch("LM_PX1",&fPositionX1,"PX1/D");
	outTree->Branch("LM_PX2",&fPositionX2,"PX2/D");
	outTree->Branch("LM_PY1",&fPositionY1,"PY1/D");
	outTree->Branch("LM_PY2",&fPositionY2,"PY2/D");
	//outTree->Branch("baseline",baseline,"baseline[1024]/F");
	outTree->Branch("time",&timeline,"time/D");
	outTree->Branch("LM_P2_Integral",&integralP2,"IP2/D");
	calibTree->Branch("LM_P2_Integral_mean",&integralP2_mean,"IP2_mean/D");
	calibTree->Branch("LM_P2_Integral_mean_error",&integralP2_mean_error,"IP2_mean_error/D");
	calibTree->Branch("LM_P2_Integral_sigma",&integralP2_sigma,"IP2_sigma/D");
	calibTree->Branch("LM_P2_Integral_sigma_error",&integralP2_sigma_error,"IP2_sigma_error/D");

    /**************************************
     * read entries
     **************************************
    */
	for (int j = 0; j < entries; ++j){
		gSystem->Sleep (sleep);

		t->GetEntry(j);

		//TRIGGER SELECTION
		if(caen.trace[flagChannel][400]>1000 && caen.trace[flagChannel][800]<3000){
			timeline = caen.system_clock;

			/**************************************
		     * Peaks estimation
		     **************************************
		    */
			for (int i = 0; i < nbins; ++i){
				g->SetPoint(i, i, caen.trace[channel][i]);
			}

			Double_t y_max = TMath::MaxElement(g->GetN(),g->GetY());
			Float_t * source = new Float_t[nbins];
			Float_t * dest = new Float_t[nbins];

			for (int i = 0; i < nbins; ++i){
				source[i]=y_max-caen.trace[channel][i];
				g->SetPoint(i, i, source[i]);
			}

			//Use TSpectrum to find the peak candidates
			TSpectrum *s = new TSpectrum();
		   	Int_t nfound = s->SearchHighRes(source, dest, nbins, 3, 2, kTRUE, 2, kFALSE, 5);

		    /**************************************
		     * Background estimation
		     **************************************
		    */
		    Int_t  	ssize = nbins;
		    Int_t  	numberIterations = 20;
		    Int_t  	direction = s->kBackIncreasingWindow;
		    Int_t  	filterOrder = s->kBackOrder2;
		    bool  	smoothing = kFALSE;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  //Hxswg::utils::multiplyGraph(   ExpLimit, XSecMELA);
  //Hxswg::utils::multiplyGraph( ExpLimitp1, XSecMELA);
  //Hxswg::utils::multiplyGraph( ExpLimitp2, XSecMELA);
 
  //Scale exclusion XSec in fb
  scaleGraph(ObsLimit  , 0.001); //pb to fb
  scaleGraph(ExpLimitm2, 0.001); //pb to fb
  scaleGraph(ExpLimitm1, 0.001); //pb to fb
  scaleGraph(ExpLimit  , 0.001); //pb to fb
  scaleGraph(ExpLimitp1, 0.001); //pb to fb
  scaleGraph(ExpLimitp2, 0.001); //pb to fb

  //scal eTH cross-section and limits according to scale factor 
  //this only apply to NarrowResonnance case
  if(strengthLimit){
     Hxswg::utils::divideGraph(ObsLimit   , THXSec);
     Hxswg::utils::divideGraph(ExpLimitm2 , THXSec);
     Hxswg::utils::divideGraph(ExpLimitm1 , THXSec);
     Hxswg::utils::divideGraph(ExpLimit   , THXSec);
     Hxswg::utils::divideGraph(ExpLimitp1 , THXSec);
     Hxswg::utils::divideGraph(ExpLimitp2 , THXSec);
     Hxswg::utils::divideGraph(THXSec     , THXSec);
  }


  //limits in terms of signal strength
  TCanvas* c = new TCanvas("c", "c",800,800);
  c->SetGridx();
  c->SetGridy();
  TH1F* framework = new TH1F("Graph","Graph",1,strengthLimit?199:199,2500); //3000);
  framework->SetStats(false);
  framework->SetTitle("");
  framework->GetXaxis()->SetTitle("M_{H} [GeV]");
  framework->GetYaxis()->SetTitleOffset(1.70);
  if(strengthLimit){
  framework->GetYaxis()->SetTitle("#mu = #sigma_{95%} / #sigma_{th}");
  framework->GetYaxis()->SetRangeUser(1E-4,1E3);
  c->SetLogy(true);
  }else{
  framework->GetYaxis()->SetTitle((string("#sigma_{95%} (") + prod +" #rightarrow H #rightarrow ZZ) (pb)").c_str());
  framework->GetYaxis()->SetRangeUser(1E-3,1E3);
  c->SetLogy(true);
  }
  framework->GetXaxis()->SetLabelOffset(0.007);
  framework->GetXaxis()->SetLabelSize(0.03);
  framework->GetXaxis()->SetTitleOffset(1.0);
  framework->GetXaxis()->SetTitleFont(42);
  framework->GetXaxis()->SetTitleSize(0.035);
  framework->GetYaxis()->SetLabelFont(42);
  framework->GetYaxis()->SetLabelOffset(0.007);
  framework->GetYaxis()->SetLabelSize(0.03);
  framework->GetYaxis()->SetTitleOffset(1.3);
  framework->GetYaxis()->SetTitleFont(42);
  framework->GetYaxis()->SetTitleSize(0.035);
  framework->Draw();

  
  TGraph* TGObsLimit   = ObsLimit;  TGObsLimit->SetLineWidth(2);
  TGraph* TGExpLimit   = ExpLimit;  TGExpLimit->SetLineWidth(2); TGExpLimit->SetLineStyle(2);
  TCutG* TGExpLimit1S  = GetErrorBand("1S", ExpLimitm1, ExpLimitp1);  
  TCutG* TGExpLimit2S  = GetErrorBand("2S", ExpLimitm2, ExpLimitp2);  TGExpLimit2S->SetFillColor(5);
  THXSec->SetLineWidth(2); THXSec->SetLineStyle(1); THXSec->SetLineColor(4);

  TGExpLimit->SetLineColor(1);  TGExpLimit->SetLineStyle(2);
  TGObsLimit->SetLineWidth(2);  TGObsLimit->SetMarkerStyle(20);
  TGExpLimit2S->Draw("fc same");

int main(int argc, char** argv){

	// because root is a piece of shit	
	gSystem->Load("libTree");

	if (argc != 2){
		printf("Usage: ./pulse_shape <gain_meas root file>\n\te.g. ./pulse_shape test.root\n");
		return -1;
	}

	char* file_name = argv[1];
//	float t_min = -70e-9, t_max = 130e-9, u_min = -1., u_max = 2.; 
	int sample_len;
	float baseline;
	

	// open root file and extract tree
 	printf("[Pulse Shape] - Opening file '%s'..\n", file_name);
	TFile* in_file = TFile::Open(file_name, "READ");
	TTree* in_tree = (TTree*)in_file->Get("outtree");
	in_tree->SetBranchAddress("len", &sample_len);
	in_tree->SetBranchAddress("baseline", &baseline);
	in_tree->GetEntry(0);

	double* u = new double[sample_len];
	double* t = new double[sample_len];
	in_tree->SetBranchAddress("amplitude", u);
	in_tree->SetBranchAddress("time", t);
	const int nEvents = in_tree->GetEntries();

 	printf("[Pulse Shape] - Found tree with %i events, %i samples per event.\n", nEvents, sample_len);

	TFile* out_file = new TFile("out.root", "RECREATE");
	TTree* out_tree = new TTree("out_tree", "outtree");
	double riseTime, fallTime, pulseDuration;
	out_tree->Branch("risetime", &riseTime);
	out_tree->Branch("falltime", &fallTime);
	out_tree->Branch("pulseduration", &pulseDuration);
	
	TCanvas* c1 = new TCanvas();
	TGraph* pulse = new TGraph();
	pulse->SetTitle("Output pulse;t [s];U [V]");
	pulse->SetMarkerStyle(7);
	TGraph* rf = new TGraph();	// drawing rise and fall time points
	rf->SetMarkerStyle(8);
	rf->SetMarkerColor(46);
	TF1* bl = new TF1("baseline", "[0]", -100, 100);	// baseline
	bl->SetLineColor(38);
	
        // loop over data
	float uMax, lPos, hPos, lTime, hTime, rTime, buf;
	int maxEntry;
        for (int iEvent = 0; iEvent < nEvents; iEvent++){
                uMax = -1.;
                in_tree->GetEntry(iEvent);
                for (int i = 0; i < sample_len; i++){
                        u[i] *= -1;
                        pulse->SetPoint(i, t[i], u[i]);
                        // find Maximum by Hand because root apparently isnt able to do so
                        if (u[i] > uMax){
                                uMax = u[i];
                                maxEntry = i;
                        }
                }
		// get 10% and 90% amplitude
                lPos = 0.1*(uMax - baseline) + baseline;
                hPos = 0.9*(uMax - baseline) + baseline;
                
		// get rise time -> start at maximum and go left
		lTime = -1;
                hTime = -1;
                for (int i = maxEntry; (buf = pulse->GetY()[i]) > 0.9*lPos; i--){
                        if ( buf > hPos )
                                hTime = pulse->GetX()[i];
                        if ( buf > lPos ){
                                lTime = pulse->GetX()[i];
                        }
                }
                riseTime = hTime - lTime;
		rf->SetPoint(0, lTime, lPos);               
		rf->SetPoint(1, hTime, hPos);               
 
		// get fall time -> start at maximum and go right
		rTime = -1;
                hTime = -1;
                for (int i = maxEntry; (buf = pulse->GetY()[i]) > 0.9*lPos; i++){
                        if ( buf > hPos )
                                hTime = pulse->GetX()[i];
                        if ( buf > lPos ){
                                rTime = pulse->GetX()[i];
                        }
                }
                fallTime = rTime - hTime;
		pulseDuration = rTime - lTime; 
		rf->SetPoint(2, rTime, lPos);               
		rf->SetPoint(3, hTime, hPos);               

		out_tree->Fill();

                // draw & save every 500th event
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    TCanvas* c3[vol_size];
    TCanvas* c4[vol_size];
    
    TMultiGraph *Cleanwaves[vol_size];
    TCanvas* expfit_longtau_c[vol_size];
    TCanvas* expfit_AP_c[vol_size];
    
    cout<<"////////////"<< endl;
    cout<<"****----->Voltage Breakdown calculation ***"<< endl;
    
   
    vector <Double_t> pe_volt;
    TGraph *Vbias_ver= new TGraph();
    
    
        //Change to not recalculate the pe
        //pe_volt.push_back(6.87435e-02);
        /*pe_volt.push_back( 1.20426e-01);
        pe_volt.push_back(1.75262e-01);
        pe_volt.push_back(2.30936e-01);
        pe_volt.push_back(2.87958e-01);*/
        //pe_volt.push_back( 3.44156e-01);
        //Double_t VBD=55.9006;
    
    //Calculate Voltage breakdown and value of pe
    for (int i=0; i<vol_size; i++) {
        pe_volt.push_back(Amplitude_calc(vol_folders.at(i).Data(), data_size));
        V_meas = vol_folders.at(i).Atof();
        Vbias_ver->SetPoint(i, pe_volt.at(i), V_meas);
    }
    
    TCanvas* ca= new TCanvas("Voltage Breakdown calculation","Voltage Breakdown calculation",100,100,900,700);
    Vbias_ver->SetTitle("Voltage Breakdown calculation");
    Vbias_ver->GetYaxis()->SetTitle("Bias Volatge [V]");
    Vbias_ver->GetYaxis()->SetTitleOffset(1.2);
    Vbias_ver->GetXaxis()->SetTitle("Mean peak amplitude [V]");
    Vbias_ver->Draw("AP*");
    ca->SetGrid();
    
    TPaveText * pv = new TPaveText(0.2,0.65,0.35,0.74,"brNDC");
    
    cout<<"////////////"<< endl;
    cout<<"****----->Voltage Breakdown fit ***"<< endl;
    
    TFitResultPtr fit = Vbias_ver->Fit("pol1","S");
    Double_t VBD= fit->Value(0);
    
    Char_t VBD_text[20];
    sprintf(VBD_text,"V_{BD} = %2.2f",VBD);
    pv->AddText(VBD_text);
    pv->Draw();
    
    if (globalArgs.save_all==1) ca->Write();
    
    cout<<"////////////"<< endl;
    cout<<"****----->Noise analysis ***"<< endl;
    cout<<"////////////"<< endl;
    
    /////////////////
    // Loop over all Voltages measured
    /////////////////
    for (int i=0; i<vol_size; i++) {
        
        //Important to reinitialize, the value color* = kOrange-11 is used to plot axis of TGraph()
        
        int color1 = kOrange-11;

void plot() {

    TGraph * FWHM = new TGraph( "plot.txt" , "%lg %lg %*lg %*lg" );
    TGraph * FWTM = new TGraph( "plot.txt" , "%lg %*lg %lg %*lg" );
    TGraph * FWFM = new TGraph( "plot.txt" , "%lg %*lg %*lg %lg" );

    TGraph * ratioTH = new TGraph();
    TGraph * ratioFH = new TGraph();
    ratioTH->SetName("ratioTH");
    ratioFH->SetName("ratioFH");

    double * x = FWHM->GetX();
    double * yH = FWHM->GetY();
    double * yT = FWTM->GetY();
    double * yF = FWFM->GetY();

    for ( int i = 0 ; i < FWHM->GetN() ; i++ ) {
        ratioTH->SetPoint( i , x[i] , yT[i]/yH[i] );
        ratioFH->SetPoint( i , x[i] , yF[i]/yH[i] );
    }
	
	ratioTH->SetTitle("Gaussianity");
    ratioTH->GetXaxis()->SetTitle("Energy Width [us]");
	ratioTH->SetMarkerStyle(20);
	ratioTH->SetMarkerColor(kViolet+2);
	ratioTH->SetMarkerSize(1.5);
	ratioTH->SetLineColor(kViolet+2);
    
	ratioFH->SetMarkerStyle(29);
	ratioFH->SetMarkerColor(kYellow+2);
	ratioFH->SetMarkerSize(1.5);
	ratioFH->SetLineColor(kYellow+2);
    ratioTH->GetYaxis()->SetRangeUser(1,6);

    FWHM->SetTitle("Resolution");
    FWHM->GetYaxis()->SetRangeUser(0,17);
    FWHM->GetYaxis()->SetTitle("[keV]");
    FWHM->GetXaxis()->SetTitle("Energy Width [us]");

    FWHM->SetMarkerStyle(33);
    FWHM->SetMarkerColor(kGreen+3);
    FWHM->SetMarkerSize(1.5);
	FWHM->SetLineColor(kGreen+3);

    FWTM->SetMarkerStyle(22);
    FWTM->SetMarkerColor(kBlue);
    FWTM->SetMarkerSize(1.5);
	FWTM->SetLineColor(kBlue);

    FWFM->SetMarkerStyle(23);
    FWFM->SetMarkerColor(kRed);
    FWFM->SetMarkerSize(1.5);
	FWFM->SetLineColor(kRed);

    TCanvas * can = new TCanvas("can","can",1);
    can->Divide(2,1);
    can->cd(1);
    gPad->SetGrid();

    FWHM->Draw("AP");
    FWTM->Draw("PSAME");
    FWFM->Draw("PSAME");
	
	TLegend * leg = new TLegend(0.24,0.77,0.42,0.88);
	leg->SetTextAlign(22);
	leg->SetTextSize(0.037);
	leg->AddEntry(FWHM,"FWHM","p");
	leg->AddEntry(FWTM,"FWTM","p");
	leg->AddEntry(FWFM,"FWFM","p");
	leg->Draw();

    can->cd(2);
	gPad->SetGrid();
    ratioTH->Draw("AP");
    ratioFH->Draw("PSAME");
  
	TLine * lineTH = new TLine( ratioTH->GetXaxis()->GetXmin() , 1.82 , ratioTH->GetXaxis()->GetXmax() , 1.82 );
	TLine * lineFH = new TLine( ratioTH->GetXaxis()->GetXmin() , 2.38 , ratioTH->GetXaxis()->GetXmax() , 2.38 );

	lineTH->SetLineWidth(2);
	lineTH->SetLineStyle(9);
	lineTH->SetLineColor(kViolet+2);

	lineFH->SetLineWidth(2);
	lineFH->SetLineStyle(9);
	lineFH->SetLineColor(kYellow+2);

	lineTH->Draw();
    lineFH->Draw();
	
	TLegend * leg2 = new TLegend(0.14,0.77,0.45,0.88);
	leg2->SetTextAlign(22);
	leg2->SetTextSize(0.037);
	leg2->AddEntry(ratioTH,"FWTM/FWHM","p");
	leg2->AddEntry(ratioFH,"FWFM/FWHM","p");
	leg2->Draw();

	return;
}

//.........这里部分代码省略.........
    else if(theory=="kidonakis"){ 
      if(variable.Contains("topPt")){
	xmin=   0.;
	xmax=1500.;
	binwidth=1.;
      }
      else if(variable.Contains("topY")){
	xmin=-3.8;
	xmax= 3.8;
	binwidth=0.01;
      }
    }
    // fill data in binned TH1F
    hist= new TH1F ( variable, variable, (int)((xmax-xmin)/binwidth), xmin, xmax);
    TH1F* ori=(TH1F*)hist->Clone("original points");
    std::cout << "fine binned theory prediction has " << hist->GetNbinsX() << " bins" << std::endl;
    std::cout << "loaded values from .dat file: " << std::endl;
    // list all data values loaded and the corresponding bin
    for(int bin=1; bin<=Nbins; ++bin){
      double x=-999;
      double y=-999;
      // NB: choose if loaded data is interpreted as points with nothing 
      //     between (like kidonakis) or as integrated over the bin range (like ahrens)
      if(points){
	// check if you are still inside the array
	//std::cout << "data point " << bin-1 << "/" << sizeof(Xvalues)/sizeof(double) << std::endl;
	if(rawHist->GetPoint(bin-1, x, y)!=-1){
	  //x=Xvalues[bin-1]; // get value from data points
	  x+=0.5*binwidth;  // add half the binwidth to get the center of the bin
	}
      }
      else{
	x=hist->GetBinCenter(bin);  // get bin center
	y=rawHist->GetY()[bin-1];
      }
      if(x!=-999){
	std::cout << "data point: " << bin;
	std::cout << "(<x>=" << x << ")-> bin";
	// get bin in target (fine binned) plot
	int bin2=bin;
	if(points) bin2=hist->FindBin(x);
	//double y=Yvalues[bin2-1];
	std::cout << bin2 << " (";
	std::cout << hist->GetBinLowEdge(bin2) << ".." << hist->GetBinLowEdge(bin2+1);
	std::cout << "): " << y << std::endl;
	// fill target (fine binned) plot
	if(!points) hist->SetBinContent(bin2, y);
	// mark bins coming from the original prediction
	ori->SetBinContent(bin2, 1.);
	// -------------------------------
	// fit range without data entries
	// -------------------------------
	// NB: needed if loaded data is interpreted as points with nothing 
	//     between (like kidonakis)
	if(points){
	  // perform a linear fit wrt previous point
	  // get the two points (this bin and the previous one)
	  int binPrev= (bin==1&&variable=="topPt") ? 0 : hist->FindBin(Xvalues[bin-2]+0.5*binwidth);
	  double x2=-1;
	  double x1= 0;
	  double y2=-1;
	  double y1= 0.;
	  rawHist->GetPoint(bin-1, x2, y2);
	  x2+=0.5*binwidth;
	  if(bin==1&&variable=="topPt"){
	    y1=0;

/*
root -l 'sigmapeak.C+(0.025)'
// 
Processing sigmapeak.C+(0.025)...
background only: bkg_nsigma_0_99 = -0.352803
sig_nsigma_0_59  = -0.608621
sig_nsigma_0_99  = 2.1472
sig_nsigma_60_99 = 4.14042
sig_nsigma_70_90 = 5.57689
sig_nsigma_75_85 = 8.056
*/
void sigmapeak(Double_t signal_area=0.025)
{
   Double_t bkg_mean = 0;
   Double_t bkg_sigma = 0.001;

   Double_t x[100];
   Double_t y[100];
   Int_t np = 100;
   
   TRandom3 rand;
   
   for (int i=0; i<np; ++i) {
      x[i] = i;
      y[i] = rand.Gaus(bkg_mean,bkg_sigma);
   }
   
   TGraph* gr = new TGraph(np,x,y);
   gr->SetNameTitle("gr",Form("#sigma = %0.1f",bkg_sigma));
   gr->SetMarkerStyle(7);
   gr->SetMarkerColor(46);
   new TCanvas();
   gr->Draw("ap");

   Double_t bkg_nsigma_0_99 = Nsigma(gr->GetY(), 0,99, bkg_sigma);
   cout<< "background only: bkg_nsigma_0_99 = " << bkg_nsigma_0_99 <<endl;

   // add signal

   Double_t signal_mean = 80;
   Double_t signal_sigma = 3;

   for (int i=0; i<gr->GetN(); ++i) {
      Double_t xx = (gr->GetX()[i] - signal_mean)/signal_sigma;
      Double_t arg = 0.5*xx*xx;
      Double_t exp = arg < 50? TMath::Exp(-arg): 0;
      Double_t signal = signal_area/(TMath::Sqrt(TMath::TwoPi())*signal_sigma) * exp;
      gr->SetPoint(i, gr->GetX()[i], gr->GetY()[i] + signal);
   }

   gr->SetTitle(Form("#sigma_{bkg} = %0.3f signal: area = %0.3f mean = %0.0f sigma = %0.1f", bkg_sigma,signal_area,signal_mean,signal_sigma));

   gr->Draw("apl");
   gr->Fit("gaus", "R", "", signal_mean - 5*signal_sigma, signal_mean+5*signal_sigma);
   Double_t fit_area = 2.5 * gr->GetFunction("gaus")->GetParameter("Constant") * gr->GetFunction("gaus")->GetParameter("Sigma");
   cout<< "Area under fitted gaussian = " << fit_area <<endl;

   // titmax();
   gPad->Modified();    // to create box (NB: the pad was not drawn yet at this point!)
   gPad->Update();
   TPaveText* tit = (TPaveText*)gPad->GetPrimitive("title");
   tit->SetX1NDC(0.);
   tit->SetX2NDC(1.);
   tit->SetY1NDC(0.9);
   tit->SetY2NDC(1.);
   gPad->Modified();    // to update the pad
   gPad->Update();

   Double_t sig_nsigma_0_59 = Nsigma(gr->GetY(), 0,59, bkg_sigma);
   cout<< "sig_nsigma_0_59  = " << sig_nsigma_0_59 <<endl;

   Double_t sig_nsigma_0_99 = Nsigma(gr->GetY(), 0,99, bkg_sigma);
   cout<< "sig_nsigma_0_99  = " << sig_nsigma_0_99 <<endl;

   Double_t sig_nsigma_60_99 = Nsigma(gr->GetY(), 60,99, bkg_sigma);
   cout<< "sig_nsigma_60_99 = " << sig_nsigma_60_99 <<endl;

   Double_t sig_nsigma_70_90 = Nsigma(gr->GetY(), 70,90, bkg_sigma);
   cout<< "sig_nsigma_70_90 = " << sig_nsigma_70_90 <<endl;

   Double_t sig_nsigma_75_85 = Nsigma(gr->GetY(), 75,85, bkg_sigma);
   cout<< "sig_nsigma_75_85 = " << sig_nsigma_75_85 <<endl;

   Double_t ys5[100];
   smooth5(np, gr->GetY(), ys5);
   TGraph* gr5 = new TGraph(np, x, ys5);
   gr5->SetNameTitle("gr5","smoothed on 5 points");
   gr5->SetMarkerStyle(7);
   gr5->SetMarkerColor(2);
   gr5->SetLineColor(2);
   new TCanvas;
   gr5->Draw("apl");

   Double_t ys7[100];
   smooth7(np, gr->GetY(), ys7);
   TGraph* gr7 = new TGraph(np, x, ys7);
   gr7->SetNameTitle("gr7","smoothed on 7 points");
   gr7->SetMarkerStyle(7);
   gr7->SetMarkerColor(8);
   gr7->SetLineColor(8);
//.........这里部分代码省略.........