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

void
plotMaxLikelihood(const char* filename, const char* channel, double min_=0., double max_=-1., bool log_=false, std::string dataset_="CMS Preliminary,  H#rightarrow#tau#tau,  4.9 fb^{-1} at 7 TeV, 19.4 fb^{-1} at 8 TeV", std::string xaxis_="m_{H} [GeV]", std::string yaxis_="best fit for #sigma/#sigma_{SM}", bool mssm_=false)
{
  TFile* file = TFile::Open(filename);
  // retrieve TGraphs from file
  TGraph* expected = (TGraph*)file->Get(std::string(channel).append("/expected").c_str());
  TGraphAsymmErrors* innerBand = (TGraphAsymmErrors*)file->Get(std::string(channel).append("/innerBand").c_str());
  // set up styles
  SetStyle();
  // create the unit line
  TGraph* unit = 0;
  if(!mssm_){
    unit = new TGraph();
    for(int imass=0, ipoint=0; imass<expected->GetN(); ++imass){
      unit->SetPoint(ipoint, expected->GetX()[imass], 1.); ++ipoint;
    }
  }
  // set proper maximum
  float max = maximum(expected, max_);
  // do the plotting 
  TCanvas canv = TCanvas("canv", "Limits", 600, 600);
  plottingLimit(canv, innerBand, 0, expected, 0, unit, xaxis_, yaxis_, min_, max, log_, "BESTFIT", "", mssm_);
  // setup CMS Preliminary
  CMSPrelim(dataset_.c_str(), "", 0.145, 0.835);
  // write results to files
  canv.Print(std::string(channel).append("_").append("bestfit").append(".png").c_str());
  canv.Print(std::string(channel).append("_").append("bestfit").append(".pdf").c_str()); 
  canv.Print(std::string(channel).append("_").append("bestfit").append(".eps").c_str()); 
  return;
}

void getStartFinFrom(const TGraph& gr,
                     Double_t& start, Double_t& fin,
                     const Bool_t filterZero=kTRUE) {
   // get a nice time axis for the domain of the graph
   start = fin = 0;
   const Int_t n = gr.GetN();
   Bool_t first = kTRUE;
   const Double_t* y = gr.GetX();
   for (Int_t i=0; i<n; ++i, ++y) {
      if (filterZero && ((*y)==0)) {
         continue;
      } else {
         if (first) {
            start = fin = *y;
            first = kFALSE;
         } else {
            if ( (*y) < start ) {
               start = *y;
            }
            if ( (*y) > fin ) {
               fin = *y;
            }
         }
      }
   }
}

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 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 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");
}

//______________________________________________________________________________
void DrawEvolution(const char* file, bool normalized=true)
{

  TFile* f = TFile::Open(gSystem->ExpandPathName(file));

  if (!f) return;

  TCanvas* c = new TCanvas("mch-status-evolution","mch-status-evolution");

  c->SetGridy();
  c->SetTicky();

  c->Draw();

  TLegend* l = new TLegend(0.1,0.7,0.3,0.95,"ch evolution");

  TGraph* g = static_cast<TGraph*>(f->Get("nbad"));
  if (!g) return;

  int runmin = TMath::Nint(g->GetX()[0]);
  int runmax = TMath::Nint(g->GetX()[g->GetN()-1]);

  cout << Form("Run range found in file %s = %d - %d",file,runmin,runmax) << endl;

  double ymax(0.4);

  TH2* h = new TH2F("hframe","hframe;Run number;Fraction of dead channels",100,runmin-200,runmax+200,100,0,ymax);

  gStyle->SetOptStat(kFALSE);
  h->Draw();
  h->GetXaxis()->SetNoExponent();
  h->GetXaxis()->SetNdivisions(505);

  gStyle->SetOptTitle(kFALSE);

  DrawPeriod(runmin,runmax,115881,117222,0,ymax,"10b");
  DrawPeriod(runmin,runmax,119159,120824,0,ymax,"10c");
  DrawPeriod(runmin,runmax,122374,126424,0,ymax,"10d");
  DrawPeriod(runmin,runmax,127724,130850,0,ymax,"10e");
  DrawPeriod(runmin,runmax,133005,134929,0,ymax,"10f");
  DrawPeriod(runmin,runmax,135658,136376,0,ymax,"10g");
  DrawPeriod(runmin,runmax,137133,139513,0,ymax,"10h");
  DrawPeriod(runmin,runmax,143856,146860,0,ymax,"11a");
  DrawPeriod(runmin,runmax,148370,150702,0,ymax,"11b");
  DrawPeriod(runmin,runmax,151566,154583,0,ymax,"11c");
  DrawPeriod(runmin,runmax,158084,159606,0,ymax,"11d");
  DrawPeriod(runmin,runmax,160677,162717,0,ymax,"11e");
  DrawPeriod(runmin,runmax,162933,165744,0,ymax,"11f");
  DrawPeriod(runmin,runmax,167703,170593,0,ymax,"11h");

  // 2012

  DrawPeriod(runmin,runmax,176661,177295,0,ymax,"12a");
  DrawPeriod(runmin,runmax,177384,178053,0,ymax,"12b");
  DrawPeriod(runmin,runmax,179603,180569,0,ymax,"12c");
  DrawPeriod(runmin,runmax,183913,186320,0,ymax,"12d");
  DrawPeriod(runmin,runmax,186365,186602,0,ymax,"12e");
  DrawPeriod(runmin,runmax,186668,188123,0,ymax,"12f");
  DrawPeriod(runmin,runmax,188362,188503,0,ymax,"12g");
  DrawPeriod(runmin,runmax,189122,190110,0,ymax,"12h");

  // 2013

  DrawPeriod(runmin,runmax,195344,195483,0,ymax,"13b");
  DrawPeriod(runmin,runmax,195529,195677,0,ymax,"13c");
  DrawPeriod(runmin,runmax,195681,195873,0,ymax,"13d");
  DrawPeriod(runmin,runmax,195949,196311,0,ymax,"13e");
  DrawPeriod(runmin,runmax,196433,197388,0,ymax,"13f");

  // 2015
  // periods are from the logbook, taking only PHYSICS% partitions
  // into account
  DrawPeriod(runmin,runmax,213329,215151,0,ymax,"15a");
  DrawPeriod(runmin,runmax,215160,215579,0,ymax,"15b");
  DrawPeriod(runmin,runmax,215580,219968,0,ymax,"15c");
  DrawPeriod(runmin,runmax,220050,223195,0,ymax,"15d");
  DrawPeriod(runmin,runmax,223227,224779,0,ymax,"15e");
  DrawPeriod(runmin,runmax,224826,226606,0,ymax,"15f");
  DrawPeriod(runmin,runmax,228936,231321,0,ymax,"15g");
  DrawPeriod(runmin,runmax,232914,234050,0,ymax,"15h");
  DrawPeriod(runmin,runmax,235196,236866,0,ymax,"15i");
  DrawPeriod(runmin,runmax,236965,238621,0,ymax,"15j");
  DrawPeriod(runmin,runmax,238890,239025,0,ymax,"15k");
  DrawPeriod(runmin,runmax,239319,241544,0,ymax,"15l");
  DrawPeriod(runmin,runmax,244340,244628,0,ymax,"15n");
  DrawPeriod(runmin,runmax,244918,246994,0,ymax,"15o");

  // 2016
  DrawPeriod(runmin,runmax,247189,247193,0,ymax,"16a");
  DrawPeriod(runmin,runmax,247723,250585,0,ymax,"16b");
  DrawPeriod(runmin,runmax,250647,252091,0,ymax,"16c");
  
  DrawPeriod(runmin,runmax,252191,252518,0,ymax,"16d");
  DrawPeriod(runmin,runmax,252598,253609,0,ymax,"16e");
  DrawPeriod(runmin,runmax,253613,253979,0,ymax,"16f");
  DrawPeriod(runmin,runmax,253991,254332,0,ymax,"16g");
  DrawPeriod(runmin,runmax,254373,255469,0,ymax,"16h");
  DrawPeriod(runmin,runmax,255515,255650,0,ymax,"16i");
  DrawPeriod(runmin,runmax,256083,256420,0,ymax,"16j");
//.........这里部分代码省略.........

// 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 msugraThExcl(TString msugrafile = "msugra_status.txt",
		  TGraph* staulsp=0, TGraph* negmasssq=0, TGraph* noRGE=0,
		  TGraph* noEWSB=0, TGraph* tachyons=0) {
  TTree msugra;
  msugra.ReadFile(msugrafile);
  
  float m0, m12;
  int status;

  msugra.SetBranchAddress("m0",&m0);
  msugra.SetBranchAddress("m12",&m12);
  msugra.SetBranchAddress("status",&status);

  TH2D *stat1hist = new TH2D("stat1hist","stat1hist",250,0,5000,100,0,1000);
  TH2D *stat2hist = new TH2D("stat2hist","stat2hist",250,0,5000,100,0,1000);
  TH2D *stat3hist = new TH2D("stat3hist","stat3hist",250,0,5000,100,0,1000);
  TH2D *stat4hist = new TH2D("stat4hist","stat4hist",250,0,5000,100,0,1000);
  TH2D *stat5hist = new TH2D("stat5hist","stat5hist",250,0,5000,100,0,1000);
  for(int i=0; i<msugra.GetEntries(); i++) {
    msugra.GetEntry(i);
    switch(status) {
    case 5:
      stat5hist->Fill(m0-1,m12-1,1);
    case 3:
      stat3hist->Fill(m0-1,m12-1,1);
    case 4:
      stat4hist->Fill(m0-1,m12-1,1);
    case 2:
      stat2hist->Fill(m0-1,m12-1,1);
      break;
    case 1:
      stat1hist->Fill(m0-1,m12-1,1);
    default: break;
    }
  }

  smoothHisto(stat1hist,3);
  smoothHisto(stat2hist,3);
  smoothHisto(stat3hist,3);
  smoothHisto(stat4hist,3);
  smoothHisto(stat5hist,3);

  stat1hist->SetContour(1);
  stat1hist->SetContourLevel(0,0.5);

  stat2hist->SetContour(1);
  stat2hist->SetContourLevel(0,0.5);

  stat3hist->SetContour(1);
  stat3hist->SetContourLevel(0,0.5);

  stat4hist->SetContour(1);
  stat4hist->SetContourLevel(0,0.5);

  stat5hist->SetContour(1);
  stat5hist->SetContourLevel(0,0.5);

  TCanvas* c = new TCanvas;

  /// DRAW STAU LSP
  stat1hist->Draw("contlist");
  c->Update();

  if (!staulsp) staulsp = new TGraph;
  TObjArray* conts_1 = (TObjArray*) gROOT->GetListOfSpecials()->FindObject("contours");
  TList* cont_1 = (TList*) conts_1->At(0);
  TGraph* mystaulsp = (TGraph*) cont_1->At(1);
  for(int pt=0; pt<mystaulsp->GetN(); pt++) {
    staulsp->SetPoint(pt,mystaulsp->GetX()[pt],mystaulsp->GetY()[pt]);
  }
  staulsp->SetPoint(staulsp->GetN(),200,1010);
  staulsp->SetPoint(staulsp->GetN(),-10,1010);
  staulsp->SetPoint(staulsp->GetN(),-10,-10);
  //staulsp->SetFillColor(kBlue);
  staulsp->SetFillColor(kMagenta-10);
  staulsp->SetLineColor(kGray+2);
  staulsp->SetLineWidth(2);

  /// DRAW -VE MASS SQUARED
  stat2hist->Draw("contlist");
  c->Update();

  if (!negmasssq) negmasssq = new TGraph;
  TObjArray* conts_2 = (TObjArray*) gROOT->GetListOfSpecials()->FindObject("contours");
  TList* cont_2 = (TList*) conts_2->At(0);
  TGraph* mynegmasssq = (TGraph*) cont_2->At(0);
  for(int pt=0; pt<mynegmasssq->GetN(); pt++) {
    negmasssq->SetPoint(pt,mynegmasssq->GetX()[pt],mynegmasssq->GetY()[pt]);
  }
  negmasssq->SetPoint(negmasssq->GetN(),5010,1010);
  negmasssq->SetPoint(negmasssq->GetN(),5010,-10);
  negmasssq->SetPoint(negmasssq->GetN(),-10,-10);
  negmasssq->SetFillColor(kAzure+1);
  negmasssq->SetLineColor(kGray+2);
  negmasssq->SetLineWidth(2);

  /// DRAW NO RGE
  stat4hist->Draw("contlist");
  c->Update();

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

void makeMuVMassPlot(bool iUseWWType = false) { 
  SetStyle();
  TCanvas *lCan = new TCanvas("A","A",600,600);
  // lCan->SetGridx(1);
  //lCan->SetGridy(1);
  lCan->SetRightMargin(0.14);

  double *lTX1 = new double[2];
  double *lTX2 = new double[2];
  double lMinNLL = 1000;
  double lVMin = 0;
  double *lMin = new double[36];     
  if(!iUseWWType) for(int i0 = 0; i0 < 36; i0++) { lMin[i0] = getMinNLL(110+i0*1.); if(lMin[i0] < lVMin) {lVMin = lMin[i0]; lTX1[0] = 110+i0*1.;}}
  //lMin[17] = (lMin[16]+lMin[18])/2.;
  //lMin[21] = (lMin[20]+lMin[22])/2.;
  //lMin[29] = (lMin[28]+lMin[30])/2.;
  //lMin[34] = (lMin[33]+lMin[35])/2.;

  TFile *lFile = new TFile("/afs/cern.ch/user/p/pharris/public/massscan/cmb+.root");
  TTree *lTree = lFile->FindObjectAny("limit");
  TH2F *lH = new TH2F("2D","2D",36,109.5,145.5,50,-2.,3.);
  float  lNLL  = 0; lTree->SetBranchAddress("nll"     ,&lNLL);
  float  lMNLL = 0; lTree->SetBranchAddress("deltaNLL",&lMNLL);
  double lMH   = 0; lTree->SetBranchAddress("mh"      ,&lMH);
  float  lR    = 0; lTree->SetBranchAddress("r"       ,&lR);
  
  if(iUseWWType) { 
    for(int i0 = 0; i0 < lTree->GetEntries(); i0++) { 
      lTree->GetEntry(i0);
      if(lR < 0.1 && lR > 0)  lMin[int(lMH-110)] = -lMNLL;
    }
    lVMin = 10000;
    for(int i0 = 0; i0 < lTree->GetEntries(); i0++) { 
      lTree->GetEntry(i0);
      double pMin =  lMin[int(lMH-110)] + lMNLL;
      if(pMin < lVMin) lVMin = pMin;
    }
  }
  for(int i0 = 0; i0 < lTree->GetEntries(); i0++) { 
    lTree->GetEntry(i0);
    //if(lMH == 125) continue;
    lNLL = 0.; //lMin  = 0.;
    lH->SetBinContent(lH->GetXaxis()->FindBin(lMH),lH->GetYaxis()->FindBin(lR),(lMNLL+lMin[lH->GetXaxis()->FindBin(lMH)-1]-lVMin)); 
    if(lMH == lTX1[0] && lMNLL < lMinNLL) {lMinNLL = lMNLL; lTX2[0] = lR;}
  }
  TH2F* lHC = lH->Clone("2D_v2");
  double lCont[3];
  lCont[0] = 1.17; 
  lCont[1] = 3.0;
  lCont[2] = 9.0;
  lHC->SetContour(2,lCont);
  //lCan->SetLogz();
  lHC->Draw("cont z list");
  lCan->Update();
  lHC->Draw("colz");

  TObjArray *lContours = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
  int lTotalConts = lContours->GetSize();
  double *lTX = new double[2]; lTX[0] = 110;  lTX[1] = 145; 
  double *lTY = new double[2]; lTY[0] = -0.5; lTY[1] = 2.5;
  TGraph *lFirst = new TGraph(2,lTX,lTY); lFirst->SetLineColor(kWhite);
  lFirst->GetXaxis()->SetRangeUser(110,148);
  lFirst->Draw("al"); lFirst->SetTitle("");
  lH->GetYaxis()->SetRangeUser(-0.5,2.5);
  lFirst->GetXaxis()->SetTitle("m_{H}[GeV]");
  lFirst->GetXaxis()->SetTitleOffset(1.0);
  lFirst->GetYaxis()->SetTitle("#mu_{best-fit}");
  lFirst->GetYaxis()->SetTitleOffset(1.2);

  lH->GetXaxis()->SetTitle("m_{H}[GeV]");
  lH->GetXaxis()->SetTitleOffset(1.0);
  lH->GetYaxis()->SetTitle("#mu_{best-fit}");
  lH->GetYaxis()->SetTitleOffset(1.2);

  lTX1[1] = lTX1[0]; lTX2[1] = lTX2[1]+0.001;
  TGraph *lSecond = new TGraph(1,lTX1,lTX2); lSecond->SetMarkerStyle(34); lSecond->SetMarkerSize(3.5);
  //lSecond->Draw("p");
  
  TLegend *lL = new TLegend(0.65,0.15,0.85,0.35); lL->SetBorderSize(0); lL->SetFillColor(0); lL->SetFillStyle(0); 
  for(i0 = 0; i0 < lTotalConts; i0++){
    pContLevel = (TList*)lContours->At(lTotalConts-1.-i0);
    // Get first graph from list on curves on this level
    std::vector<double> lX; 
    std::vector<double> lY; 
    pCurv = (TGraph*)pContLevel->First();
    for(int i1 = 0; i1 < pContLevel->GetSize(); i1++){
      for(int i2  = 0; i2 < pCurv->GetN(); i2++) {lX.push_back(pCurv->GetX()[i2]); lY.push_back(pCurv->GetY()[i2]);}
      //pCurv->GetPoint(0, x0, y0);
      pCurv->SetLineColor(kBlack);//kGreen+i0);
      pCCurv = (TGraph*)pCurv->Clone();
      if(i0 == 0) pCCurv->SetFillColor(0);
      if(i0 == 1) pCCurv->SetFillColor(0);
      //if(i0 == 1) pCCurv->SetLineStyle(kDashed);
      pCCurv->SetLineWidth(3);
      pCCurv->GetXaxis()->SetRangeUser(0,3.0);
      //if(i0 == 0) pCCurv->Draw("AL");
      //if(i0 != -10) pCCurv->Draw("LC");
      //l.DrawLatex(x0,y0,val);
      pCurv = (TGraph*)pContLevel->After(pCurv); // Get Next graph
    }
//.........这里部分代码省略.........

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
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    
    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;
        int color2 = kOrange-11;
        int color3 = 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;
}

void createNNLOplot(TString theory="ahrens")
{
  // theory="ahrens", "kidonakis"
  TString outputRootFile="test.root";
  // NB: add new datset name here
  if(theory.Contains("ahrens")   ){
    outputRootFile="AhrensNLONNLL";
    //if(theory.Contains("mtt")    ) outputRootFile+="mttbar" ;
    //else if(theory.Contains("pt")) outputRootFile+="pTttbar";
    outputRootFile+=".root";
  }
  else if(theory=="kidonakis") outputRootFile="KidonakisAproxNNLO.root";
  // general options
  gStyle->SetOptStat(0);
  bool usequad=true;
  bool divideByBinwidth=true;
  // list of variables
  std::vector<TString> xSecVariables_, xSecLabel_;
  // NB: add variables for new datset name here
  TString xSecVariablesUsedAhrens[]    ={"ttbarMass", "ttbarPt"};
  TString xSecVariablesUsedKidonakis[] ={"topPt"    , "topY"   };
  if(     theory.Contains("ahrens")   ) xSecVariables_ .insert( xSecVariables_.begin(), xSecVariablesUsedAhrens   , xSecVariablesUsedAhrens    + sizeof(xSecVariablesUsedAhrens   )/sizeof(TString) );
  else if(theory=="kidonakis") xSecVariables_ .insert( xSecVariables_.begin(), xSecVariablesUsedKidonakis, xSecVariablesUsedKidonakis + sizeof(xSecVariablesUsedKidonakis)/sizeof(TString) );
  // get variable binning used for final cross sections 
  std::map<TString, std::vector<double> > bins_=makeVariableBinning();
  //std::vector<double> tempbins_;
  //double ttbarMassBins[]={345,445,545,745, 1045};
  //tempbins_.insert( tempbins_.begin(), ttbarMassBins, ttbarMassBins + sizeof(ttbarMassBins)/sizeof(double) );
  //bins_["ttbarMass"]=tempbins_;

  // loop all variables
  for(unsigned var=0; var<xSecVariables_.size(); ++var){
    TString variable=xSecVariables_[var];
    std::cout << "----------" << std::endl;
    std::cout << theory << ": " << variable << std::endl;
    // get bin boundaries
    double low =bins_[variable][0];
    double high=bins_[variable][bins_[variable].size()-1];
    // --------------------------
    // get raw nnlo theory points
    // --------------------------
    // NB: add new datset file names here
    TGraph * rawHist;
    TString predictionfolder="/afs/naf.desy.de/group/cms/scratch/tophh/CommonFiles/";
    if(theory=="ahrens"   ){
      if(variable.Contains("ttbarMass")) rawHist= new TGraph(predictionfolder+"AhrensTheory_Mttbar_8000_172.5_NLONNLL_norm.dat");//AhrensTheory_Mtt_7000_172.5_Mtt_fin.dat
      if(variable.Contains("ttbarPt"  )) rawHist= new TGraph(predictionfolder+"AhrensTheory_pTttbar_8000_172.5_NLONNLL_abs.dat");//AhrensTheory_pTttbar_7000_172.5_NLONNLL_abs.dat
    }
    else if(theory=="kidonakis"){
      if(variable.Contains("topPt")) rawHist= new TGraph(predictionfolder+"pttopnnloapprox8lhc173m.dat");//"ptnormalNNLO7lhc173m.dat" //"pttopnnloapprox8lhc173m.dat"
      if(variable.Contains("topY" )) rawHist= new TGraph(predictionfolder+"ytopnnloapprox8lhc173m.dat" );//"ynormalNNLO7lhc173m.dat"  //"ytopnnloapprox8lhc173m.dat"
    }
    // NB: say if points should be interpreted as single points with 
    //     nothing in between or as integrated value over the range
    bool points=true;
    if(theory.Contains("ahrens")) points=false;
    else if(theory=="kidonakis")  points=true;
    std::cout << "input: " << rawHist->GetTitle() << std::endl;
    std::cout << "interprete values as points?: " << points << std::endl;

    // --------------------
    // convertion to TH1F*
    // --------------------
    double *Xvalues = rawHist->GetX();
    int Nbins=rawHist->GetN();
    //double *Yvalues = rawHist->GetY(); // not working
    double xmax=-1;
    double xmin=-1;
    double binwidth=-1;
    //
    TH1F* hist;
    // NB: add additional binning for new theory here
    // NB: if loaded data should be interpreted as points with 
    //     nothing in between (like kidonakis), make suree you 
    //     choose a binning that is fine enough for the 
    //     data points loaded
    if(theory=="ahrens"){
      if(variable.Contains("ttbarMass")){
	xmin= 345.;
	xmax=2720.;
	binwidth=25.;// 5 for 8TeV, 25 for 7TeV
	if(TString(rawHist->GetTitle()).Contains("8000")) binwidth=5.;
      }
      else if(variable.Contains("ttbarPt")){
	xmin= 0.;
	xmax=1300.;
	binwidth=5.;
      }
    }
    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;
      }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    double significance=0;
    double significance_help=0;
    for(int nbin=0; nbin<h_bbA->GetNbinsX()+1; nbin++){ 
      if(h_background->GetBinContent(nbin)+h_bbA->GetBinContent(nbin)!=0) significance=h_bbA->GetBinContent(nbin)/sqrt(h_background->GetBinContent(nbin)+h_bbA->GetBinContent(nbin));
      else significance=0;
      //if(significance>significance_help){
      significance_help+=significance;
      //}
    }
    std::cout << "significance " << significance_help << std::endl;
    signi->SetPoint(idx, mA, significance_help);

    //destruct TH1Fs
    h_ggA->~TH1F();
    h_ggH->~TH1F(); 
    h_ggh->~TH1F();
    h_bbA->~TH1F();
    h_bbH->~TH1F(); 
    h_bbh->~TH1F();

    std::cout << std::endl;
    std::cout << std::endl;
  }			    

  /// do the drawing
  TCanvas* canv1 = new TCanvas("canv1", "xs*BR", 600, 600);
  canv1->cd();
  //canv1->SetGridx(1);
  //canv1->SetGridy(1);
  canv1->SetLeftMargin(0.15);
  canv1->SetLogy(1);

  // draw a frame to define the range
  TH1F* hr=canv1->DrawFrame(ggA->GetX()[0]-.01, 0.0001, ggA->GetX()[ggA->GetN()-1]+.01, 10000);
  // format x axis
  hr->SetXTitle("m_{A}");
  hr->GetXaxis()->SetLabelSize(0.045);
  hr->GetXaxis()->SetTitleFont(62);
  hr->GetXaxis()->SetLabelFont(62);
  hr->GetXaxis()->SetTitleColor(1);
  hr->GetXaxis()->SetTitleOffset(1.08);
  // format y axis
  hr->SetYTitle("#sigma*BR and acc*signal_{eff}*lumi");
  hr->GetXaxis()->SetTitleFont(62);
  hr->GetYaxis()->SetLabelFont(62);
  hr->GetYaxis()->SetLabelSize(0.045);
  hr->GetYaxis()->SetTitleOffset(1.6);
  hr->GetXaxis()->SetTitleColor(1);
  hr->SetNdivisions(505, "X");

  //ggA->Draw("PLsame");
  ggA->SetLineStyle(3.);
  ggA->SetLineWidth(2.); 
  ggA->SetLineColor(kBlue);
  ggA->SetMarkerStyle(23);
  ggA->SetMarkerSize(0);
  ggA->SetMarkerColor(kBlue);

  //ggH->Draw("PLsame");
  ggH->SetLineStyle(2.);
  ggH->SetLineWidth(2.); 
  ggH->SetLineColor(kBlue);
  ggH->SetMarkerStyle(22);
  ggH->SetMarkerSize(0);
  ggH->SetMarkerColor(kBlue);

/*
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);
//.........这里部分代码省略.........