C++ TGraphErrors::GetPoint方法代码示例

void JEC_fit_Uncertainty(int N)
{
  TF1 *func[1000];
  TFile *inf = new TFile("L3Graphs_test_Icone5.root");
  TGraphErrors *g = (TGraphErrors*)inf->Get("Correction_vs_CaloPt");
  TGraphErrors *vg[1000];
  int i,k;
  double x[20],y[20],ex[20],ey[20];
  double vx[20],vy[20],vex[20],vey[20];
  for(i=0;i<g->GetN();i++)
    {
      g->GetPoint(i,x[i],y[i]);
      ex[i]=g->GetErrorX(i);
      ey[i]=g->GetErrorY(i); 
    }  
  TRandom *rnd = new TRandom();
  rnd->SetSeed(0);
  for(k=0;k<N;k++)
    {
      for(i=0;i<g->GetN();i++)
        {	
          vx[i] = rnd->Gaus(x[i],ex[i]);
          //vx[i] = x[i];
          vy[i] = rnd->Gaus(y[i],ey[i]);
          vex[i] = ex[i];
          vey[i] = ey[i];
        }
      vg[k] = new TGraphErrors(g->GetN(),vx,vy,vex,vey);
      func[k] = new TF1("func","[0]+[1]/(pow(log10(x),[2])+[3])",1,2000);
      func[k]->SetParameters(1,3,6,5);
      vg[k]->Fit(func[k],"RQ");     	
    }
  
  TCanvas *c = new TCanvas("c","c");
  gPad->SetLogx();
  g->SetMarkerStyle(20);
  g->SetMaximum(3.5);
  g->Draw("AP");
  for(k=0;k<N;k++)
    {
      func[k]->SetLineColor(5);
      func[k]->SetLineWidth(1);
      cout<<func[k]->GetChisquare()<<endl;
      vg[k]->SetMarkerColor(2);
      vg[k]->SetLineColor(2);
      vg[k]->SetMarkerStyle(21);
      //if (func[k]->GetChisquare()<0.1)
        //vg[k]->Draw("sameP");
      func[k]->Draw("same");  
    }  	 
}  

void drawTimeDifference (TDirectory* directory, TH1* refHisto, const char* fname=0)
{
  TGraphErrors* graphX = (TGraphErrors*)directory->Get("x");
  if ( graphX==0 )  return;
  TH1I* hFirst = (TH1I*)directory->Get("firstTime");
  TH1I* hLast = (TH1I*)directory->Get("lastTime");
  if ( hFirst==0 || hLast==0 )  return;

  std::string fullName("cDeltaT");
  if ( fname )  fullName += fname;
  else  fullName += directory->GetName();
  TCanvas* c = new TCanvas(fullName.c_str(),fullName.c_str());
  TH1* h = refHisto->Clone("DeltaT");
  h->Reset();
  h->SetTitle("DeltaT");

  TGraph* graph = new TGraph();
  graph->SetName("gDeltaT");

  double xg,yg;
  for ( unsigned int i=1; i<=hFirst->GetNbinsX(); ++i ) {
    std::time_t t1 = hFirst->GetAt(i);
    std::time_t t2 = hLast->GetAt(i);
    TTimeStamp ts1(hFirst->GetAt(i));
    std::cout << "Fit started at " << ts1.AsString() << std::endl;
    graphX->GetPoint(i-1,xg,yg);
    graph->SetPoint(i-1,xg,difftime(t2,t1));
  }

  double xmin,xmax,ymin,ymax;
  graph->ComputeRange(xmin,ymin,xmax,ymax);
  h->SetMinimum(0.);
  h->SetMaximum((ymax+ymin)/2.+2.*(ymax-ymin)/2.);
  h->Draw();
  graph->SetMarkerStyle(20);
//   graph->SetMarkerColor(2);
//   graph->SetLineColor(2);
  graph->Draw("P");
}

void drawEventDifference (TDirectory* directory, TH1* refHisto, const char* fname=0)
{
  TGraphErrors* graphX = (TGraphErrors*)directory->Get("x");
  if ( graphX==0 )  return;
  TH1I* hFirst = (TH1I*)directory->Get("firstEvent");
  TH1I* hLast = (TH1I*)directory->Get("lastEvent");
  if ( hFirst==0 || hLast==0 )  return;

  std::string fullName("cDeltaE");
  if ( fname )  fullName += fname;
  else  fullName += directory->GetName();
  TCanvas* c = new TCanvas(fullName.c_str(),fullName.c_str());
  TH1* h = refHisto->Clone("DeltaE");
  h->Reset();
  h->SetTitle("DeltaE");

  TGraph* graph = new TGraph();
  graph->SetName("gDeltaE");

  double xg,yg;
  for ( unsigned int i=1; i<=hFirst->GetNbinsX(); ++i ) {
    int e1 = hFirst->GetAt(i);
    int e2 = hLast->GetAt(i);
    graphX->GetPoint(i-1,xg,yg);
    graph->SetPoint(i-1,xg,e2-e1);
  }

  double xmin,xmax,ymin,ymax;
  graph->ComputeRange(xmin,ymin,xmax,ymax);
  h->SetMinimum(0.);
  h->SetMaximum((ymax+ymin)/2.+2.*(ymax-ymin)/2.);
  h->Draw();
  graph->SetMarkerStyle(20);
//   graph->SetMarkerColor(2);
//   graph->SetLineColor(2);
  graph->Draw("P");
}

//.........这里部分代码省略.........
	  MC->Sumw2();	  
	  
	  std::vector<TH2 *> systMC;
	  for(size_t isyst=0; isyst<nSysts; isyst++)
	    {
	      TH2F *h=(TH2F *)inSystF->Get(systList[isyst]+"/"+ch[ich]+"_"+profs[i]) ;
	      h->Scale(totalTTbar/h->Integral());
	      h->Add(bckgMC);
	      h->Sumw2();
	      h->SetDirectory(0);
	      systMC.push_back(h);
	    }
		

	  //get the data
	  TH2 *Data = (TH2 *) inF->Get("data/"+ch[ich]+"_"+profs[i]);            
	  Data->SetDirectory(0);
	  Data->Sumw2();

	  TGraphErrors *MCProf   = new TGraphErrors(MC->ProfileX());    MCProf->SetMarkerStyle(24);   MCProf->SetFillStyle(0); MCProf->SetName(ch[ich]+profs[i]+"mc");
	  TGraphErrors *DataProf = new TGraphErrors(Data->ProfileX());  DataProf->SetMarkerStyle(20); DataProf->SetFillStyle(0); DataProf->SetName(ch[ich]+profs[i]+"data");
	  
	  //build data/MC scale factors
	  std::vector<TGraphErrors *> data2mcProfs;
	  for(size_t isyst=0; isyst<=systMC.size(); isyst++)
	    {
	      TGraphErrors *prof= (isyst==0 ? MCProf : new TGraphErrors(systMC[isyst-1]->ProfileX()));
	      TString baseName(ch[ich]+profs[i]);
	      if(isyst) baseName += systList[isyst];
	      prof = (TGraphErrors *) prof->Clone(baseName+"data2mc");
	      for(int ip=0; ip<DataProf->GetN(); ip++)
		{
		  Double_t x,y,ydata,y_err,ydata_err;
		  prof->GetPoint(ip,x,y);         y_err=prof->GetErrorY(ip);
		  DataProf->GetPoint(ip,x,ydata); ydata_err=DataProf->GetErrorY(ip);
		  if(y<=0) continue;
		  prof->SetPoint(ip,x,ydata/y);
		  prof->SetPointError(ip,0,sqrt(pow(ydata*y_err,2)+pow(ydata_err*y,2))/pow(y,2));
		}
	      prof->SetFillColor(systColors[isyst]);
	      prof->SetFillStyle(3001+isyst%2);
	      prof->SetTitle( isyst==0 ? nomTTbarTitle : systLabels[isyst-1] );
	      //prof->SetMarkerStyle(24);   prof->SetFillStyle(3001); prof->SetMarkerColor(1+isyst); prof->SetMarkerColor(1+isyst); prof->SetLineColor(1+isyst);
	      data2mcProfs.push_back(prof);
	    }

	  
	  TH1D *MCProjY=MC->ProjectionY();
	  MCProjY->Scale(1./MCProjY->Integral()); 
	  TGraphErrors *MCProj   = new TGraphErrors(MCProjY);    MCProj->SetMarkerStyle(24); MCProj->SetFillStyle(0);  MCProj->SetName(ch[ich]+profs[i]+"projmc");
	  TH1D *DataProjY=Data->ProjectionY();
	  DataProjY->Scale(1./DataProjY->Integral());
	  TGraphErrors *DataProj = new TGraphErrors(DataProjY);  DataProj->SetMarkerStyle(20); DataProj->SetFillStyle(0); DataProj->SetName(ch[ich]+profs[i]+"projdata");
	  MCProj->SetLineColor(colors[i]);   MCProj->SetMarkerColor(colors[i]);   MCProj->SetFillColor(colors[i]); MCProj->SetFillStyle(1001);
	  DataProj->SetLineColor(colors[i]); DataProj->SetMarkerColor(colors[i]);   DataProj->SetFillColor(colors[i]);
	  
	  TPad *p=(TPad *)cproj->cd();
	  p->SetLeftMargin(0.15);
	  p->SetRightMargin(0.02);
	  p->SetTopMargin(0.05);
	  p->SetLogy();
	  MCProj->SetFillStyle(0);
	  MCProj->Draw(i==0 ? "al3" : "l3");
	  MCProj->GetYaxis()->SetRangeUser(1e-5,1.0);
	  MCProj->GetXaxis()->SetTitle( MC->GetYaxis()->GetTitle() );
	  MCProj->GetYaxis()->SetTitle( "1/N dN/dx" );

//.........这里部分代码省略.........
	while (ievt < startEv) {
		fread((void *) &event_data, 1, sizeof(event_data), fdata);
		if (feof(fdata))
			break;
		ievt++;
	}
	// filling
	ievt = 1;
	Int_t flagEnd = 0;
	Double_t chtmp;
	Double_t PedVal, itmp, Ch0Val;
	// loop on events
	cout << endl << " --- read DATA file:" << fdata << endl;
	while (ievt <= nevt && !flagEnd) {
		fread((void *) &event_data, 1, sizeof(event_data), fdata);
		if (feof(fdata))
			flagEnd = 1;
		if (ievt % (nevt / 10 + 1) == 0)
			cout << "*" << endl;
		p = (struct channel_struct *) &event_data.ch[0]; // read bunch of data
		dep = (struct channel_struct *) &event_data.ch[1]; // read bunch of data

		TGraphErrors *grCh0 = new TGraphErrors(DOMINO_NCELL);

		// loop on channels
		for (int h = 0; h < DOMINO_NCH; h++) {
			// loop on cells
			distch = (TH1F *) DistChList->At(h);
			distchsub = (TH1F *) DistChSubList->At(h);
			grPed = (TGraphErrors *) grPedList->At(h);
			distch0sub = (TH1F *) DistCh0SubList->At(h);
			if(h==0) {
				for(i = 0; i < DOMINO_NCELL;i++) {
					grPed->GetPoint(i, itmp, PedVal);
					chtmp = (Double_t)(p->data[i]);
					chtmp = chtmp - PedVal;
					grCh0->SetPoint(i,itmp, chtmp);
				}
			}
			for (int i = 0; i < DOMINO_NCELL; i++) {
				// Read pedestal value for this cell
				grPed->GetPoint(i, itmp, PedVal);
				grCh0->GetPoint(i, itmp, Ch0Val);
				//                cout << itmp << ", " << PedVal << endl;
				// Read calibration correction for this cell
				//                CalFact =

				//charge distribution for each cell, pedestal subtracted
				chtmp = (Double_t)(p->data[i]); // data value
				//                cout  << "tcell, tcell, depth: " << chtmp << ","  << p->data[i] << "," << deptmp << endl;
				distch->Fill(chtmp);
				// Check data value: must be within DOMINO Depth
				//                if(chtmp > DominoDepthADC)
				//                    cout << " === WARNING!!! Channel " << h << " Cell " << i << " has value " << chtmp << endl;
				//                cout << "Charge: " << p->data[i] << endl;
				((TH1 *) hCellList->At(h * DOMINO_NCELL + i))->Fill(chtmp);
				// Now the pedestal is subtracted
				chtmp = chtmp - PedVal;
				distchsub->Fill(chtmp);
				((TH1 *) hCellSubList->At(h * DOMINO_NCELL + i))->Fill(chtmp);
				chtmp = chtmp - Ch0Val;
				distch0sub->Fill(chtmp);
			}
			p++; // next channel
		}
		ievt++; // next event

int plotResoVsIC(){

  SetTdrStyle();


  const unsigned nIC = 10;
  const unsigned ICval[nIC] = {0,1,2,3,4,5,10,15,20,50};


  std::ostringstream label;
  TFile *fcalib[nIC];
  
  TGraphErrors *constant = new TGraphErrors();
  constant->SetName("constant");
  constant->SetTitle(";intercalib. smearing");
  constant->SetMarkerStyle(20);
  constant->SetMarkerColor(1);
  constant->SetLineColor(1);
  TGraphErrors *constantSR7 =  new TGraphErrors();
  constantSR7->SetName("constantSR7");
  constantSR7->SetTitle(";intercalib. smearing");
  constantSR7->SetMarkerStyle(23);
  constantSR7->SetMarkerColor(2);
  constantSR7->SetLineColor(2);

  TGraphErrors *noise = (TGraphErrors *) constant->Clone("noise");
  TGraphErrors *sampling = (TGraphErrors *) constant->Clone("sampling");
  TGraphErrors *samplingSR7 = (TGraphErrors *) constantSR7->Clone("samplingSR7");

  TCanvas *mycReso = new TCanvas("mycReso","mycReso",1500,1000);
  mycReso->Divide(2,5);
  TCanvas *mycR = new TCanvas("mycR","Sampling",1500,1000);
  TCanvas *mycC = new TCanvas("mycC","Constant",1500,1000);
  TCanvas *mycN = new TCanvas("mycN","Noise",1500,1000);

  gStyle->SetOptFit(1111);
  gStyle->SetOptStat(0);

  gStyle->SetStatW(0.2);
  gStyle->SetStatH(0.5);

  TLatex lat;
  char buf[500];

  TGraphErrors *gr[nIC][2];
  double x0,y0;
  double x0_7,y0_7;

  for (unsigned ic(0);ic<nIC;++ic){//loop on intercalib
    label.str("");
    label << "PLOTS/CalibReso";
    label << "_vsE";
    label << "_IC" << ICval[ic];
    label << ".root";
    fcalib[ic] = TFile::Open(label.str().c_str());
    if (!fcalib[ic]) {
      std::cout << " -- failed to open file: " << label.str() << std::endl;
      continue;
    }
    else {
      std::cout << " -- file " << label.str() << " successfully opened." << std::endl;
    }
    fcalib[ic]->cd("SR2");
    gr[ic][0] = (TGraphErrors *)gDirectory->Get("resoRecoFit2eta21pu1");
    fcalib[ic]->cd("SR7");
    gr[ic][1] = (TGraphErrors *)gDirectory->Get("resoRecoFit7eta21pu1");


    TF1 *fit = gr[ic][0]->GetFunction("reso");
    TF1 *fit7 = gr[ic][1]->GetFunction("reso");
    mycReso->cd(ic+1);
    gr[ic][0]->Draw("APE");
    fit->SetLineColor(6);
    fit->Draw("same");
    lat.SetTextSize(0.1);
    sprintf(buf,"Single #gamma, #eta=2.1, 3#times3 cm^{2}");
    lat.DrawLatexNDC(0.2,0.8,buf);
    sprintf(buf,"ICsmear = %d %%",ICval[ic]);
    lat.DrawLatexNDC(0.2,0.7,buf);

    double cval = sqrt(pow(fit->GetParameter(1),2)-pow(y0,2));
    constant->SetPoint(ic,ICval[ic]/100.,cval);
    constant->SetPointError(ic,0,fit->GetParameter(1)*fit->GetParError(1)/cval);
    noise->SetPoint(ic,ICval[ic]/100.,fit->GetParameter(2));
    noise->SetPointError(ic,0,fit->GetParError(2));
    sampling->SetPoint(ic,ICval[ic]/100.,fit->GetParameter(0));
    sampling->SetPointError(ic,0,fit->GetParError(0));
    cval = sqrt(pow(fit7->GetParameter(1),2)-pow(y0_7,2));
    constantSR7->SetPoint(ic,ICval[ic]/100.,cval);
    constantSR7->SetPointError(ic,0,fit7->GetParameter(1)*fit7->GetParError(1)/cval);    
    //constantSR7->SetPoint(ic,ICval[ic]/100.,fit7->GetParameter(1));
    //constantSR7->SetPointError(ic,0,fit7->GetParError(1));
    samplingSR7->SetPoint(ic,ICval[ic]/100.,fit7->GetParameter(0));
    samplingSR7->SetPointError(ic,0,fit7->GetParError(0));

    if (ic==0) {
      constant->GetPoint(0,x0,y0);
      constantSR7->GetPoint(0,x0_7,y0_7);
      cval = sqrt(pow(fit->GetParameter(1),2)-pow(y0,2));
      constant->SetPoint(ic,ICval[ic]/100.,cval);
//.........这里部分代码省略.........

double FitCurve(TGraphErrors* g, int debug=0)
{
  
  double vdep=0;
  if(!g) return vdep;
  
  
  TF1* f3 = new TF1("fp1", "pol1", 20, 360);
  f3->SetLineColor(4);
  
  double y;
  double xlow=350;
  double ymax=0;
  double vdep1=0;
  // prendre le point le plus haut est robuste et marche bien pour les hauts vdep
  TGraphErrors *gsmooth = MedianFilter(g);
  gsmooth = MedianFilter(gsmooth);
  for(int ipt=0; ipt<gsmooth->GetN(); ipt++)
  {
    gsmooth->GetPoint(ipt, xlow, y);
	if(y>ymax) {ymax=y; vdep1=xlow;}
  }
  
  //g=gsmooth;
  int npt = g->GetN()-5;
  xlow=350;
  double chi2=0;
  int status = 0;
  while(chi2<5. && xlow>100 && npt>=0)
  {
    g->GetPoint(npt, xlow, y);
    f3->SetRange(xlow, 350);
    status = g->Fit("fp1", "rqn");
    chi2 = f3->GetChisquare()/f3->GetNDF();
    if(debug>=1) cout<<"xlow "<<xlow<<" chi2 "<<chi2<<endl;
	npt--;
  }
  g->GetPoint(npt+2, xlow, y);
  f3->SetRange(xlow, 350);
  g->Fit("fp1", "rqn");
  vdep = xlow;
  
  //if(vdep>230) vdep = vdep1;
  //else if(fabs(vdep-vdep1)>40) vdep = vdep1;
  //vdep=vdep1;

  cout<<" Vdepl = "<<vdep<<endl;
  
  
  
  TF1* f1 = new TF1("fitsigmo", fitsigmo, 50, 350, 5);
  f1->SetParameter(0, 3.);
  f1->SetParameter(1, 0.02);
  f1->SetParameter(2, vdep);
  f1->SetParameter(3, -10.);
  f1->SetParameter(4, -0.00001);
  f1->SetParLimits(2, 100., 300.);
  f1->SetParLimits(4, -0.1, 0);
  
  
  
  
  TF1* f2 = new TF1("fitpol", fitpol, 20, 360, 6);
  f2->SetParLimits(4, 100., 300.);
  f2->SetParLimits(5, -0.1, 0);
  
 /* xlow = 30;
  chi2 = 100;
  double xmin = xlow;
  double chi2min = chi2;
  string opt="rq";
  if(debug>=2) opt="r";
  while(chi2>0.1 && xlow<vdep-80)
  {
    f2->SetParameter(0, 1.);
    f2->SetParameter(1, 0.005);
    f2->SetParameter(2, 0.0001);
    f2->SetParameter(3, 0.);
    f2->SetParameter(4, vdep);
    f2->SetParameter(5, -0.00001);
    f2->SetRange(xlow, 350);
    status = g->Fit("fitpol", opt.c_str());
    if(status!=0) status = g->Fit("fitpol", opt.c_str());
    chi2 = f2->GetChisquare()/f2->GetNDF();
    if(debug>=1) cout<<"xlow "<<xlow<<" chi2 "<<chi2<<endl;
    if(debug>=2) cout<<" chi2 "<<f2->GetChisquare()<<" ndf "<<f2->GetNDF()<<endl;
	if(chi2<chi2min) {chi2min=chi2; xmin=xlow;}
	xlow+=10;
  }
  f2->SetParameter(0, 1.);
  f2->SetParameter(1, 0.005);
  f2->SetParameter(2, 0.0001);
  f2->SetParameter(3, 0.);
  f2->SetParameter(4, vdep);
  f2->SetParameter(5, -0.00001);
  f2->SetRange(xmin, 350);
  status = g->Fit("fitpol", opt.c_str());
  if(status!=0) status = g->Fit("fitpol", opt.c_str());
  vdep = f2->GetParameter(4);
  cout<<" Vdepl2 = "<<vdep<<" xmin "<<xmin<<endl;
//.........这里部分代码省略.........

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 AnalyzeWaveforms(char *WaveformsFile = "Waveforms.root", const int nAddedChannels = 5) {
	
	//try to access waveforms file and in case of failure return
	if(gSystem->AccessPathName(WaveformsFile,kFileExists)) {
		cout << "Error: file " << WaveformsFile << " does not exsist. Run .x WaveformsFileMaker.C to create it" << endl;
		return;
	}
	//	gStyle->SetOptFit(111);
	//	gStyle->SetStatFormat("1.3E");
	//	gStyle->SetFitFormat("1.3E");
	
	// fetch the list of trees contained in the waveforms file
	// for every tree generate a waveform graph

	TFile *f = TFile::Open(WaveformsFile);
	TList *listOfKeys = f->GetListOfKeys();
	Int_t numberOfKeys = listOfKeys->GetEntries();
	TList *listOfGraphs = new TList();
	
	// if the waveform file name begins with the string "comparator" it goes in this list
	TList *listOfCompWaves = new TList();
	// if the waveform file name begins with the string "sum output" it goes in this list
	TList *listOfAdderWaves = new TList();

	for(Int_t i = 0; i < numberOfKeys; i++) {
		TString *keyName = new TString(listOfKeys->At(i)->GetName());
		TTree *tree = (TTree*)f->Get(keyName->Data());
		Float_t x = 0;
		Float_t y = 0;
		tree->SetBranchAddress("x",&x);
		tree->SetBranchAddress("y",&y);
		Int_t nentries = tree->GetEntries();

		TString *gName = new TString(keyName->Data());
		gName->Append(" graph");
		TGraphErrors *gWave = new TGraphErrors(nentries);
		gWave->SetName(gName->Data());
		gWave->SetTitle(gName->Data());
		gWave->GetXaxis()->SetTitle("Time");
		gWave->GetYaxis()->SetTitle("Voltage");
	//	gWave->SetBit(TH1::kCanRebin);

		for (Int_t j = 0; j < nentries; j++) {
			tree->GetEntry(j);
			gWave->SetPoint(j,x,y);
		}

		listOfGraphs->Add(gWave);

		if(keyName->BeginsWith("comparator"))
			listOfCompWaves->Add(gWave);

		if(keyName->BeginsWith("sum output"))
			listOfAdderWaves->Add(gWave);

	/*	TString *cName = new TString(keyName->Data());
		cName->Append(" canvas");
		TCanvas *cy = new TCanvas(cName->Data(),cName->Data(),800,600);
		gWave->Draw("AL"); */
	}
	
	cout << listOfAdderWaves->GetEntries() << endl;

	// analysis for waves with no delay

	// global variables

	Double_t xMin,xMax,yStart,yEnd;
	Int_t graphPoints;
	Double_t step;

	// comparator outputs waves sum

	TGraphErrors *gFirstCompWave = (TGraphErrors *)listOfCompWaves->First();

	graphPoints = gFirstCompWave->GetN();
	gFirstCompWave->GetPoint(0,xMin,yStart);
	gFirstCompWave->GetPoint(graphPoints - 1,xMax,yEnd); 

	step = (xMax - xMin)/graphPoints;
	cout << gFirstCompWave->GetName() << endl;
	cout << "xMin = " << xMin << "  xMax = " << xMax << "  graphPoints = " << graphPoints << endl;

	TGraphErrors *gCompSum = new TGraphErrors(graphPoints);
	gCompSum->SetLineColor(kBlue);
	gCompSum->SetLineWidth(2);
	gCompSum->SetName("Comparator Outputs Sum");
	gCompSum->SetTitle("Comparator Outputs Sum");
	Int_t nCompWaves = listOfCompWaves->GetEntries();
	Float_t gx,gy = 0;
	
	// Alpha coefficiens are now written "hard coded" here
	Float_t alphaArray[3] = {0.199,0.201,0.197};
	
	// Deleays coming from the multiplexer are written "hard coded" here
	Float_t muxDelayArray[3] = {0,77.14E-12,192.01E-12};

	for(Int_t i = 0; i < graphPoints; i++) {
		for(Int_t j = 0; j < nCompWaves; j++) {
			TGraphErrors *gCompWave = (TGraphErrors *)listOfCompWaves->At(j);
//.........这里部分代码省略.........

void CalcPeriod(char *DataFile = "drs4_peds_5buffers.dat", Int_t nevt,
		Int_t startEv = 1, char *PedFile) {

	// create progress bar
	TGHProgressBar *gProgress = ProgressBar("Calcolo periodo");

	// Redefine DOMINO Depth in ADC counts
	const Float_t DominoDepthADC = pow(2, DOMINO_DEPTH);

	// open file

	FILE *fdata = OpenDataFile(DataFile);
	struct channel_struct *p;
	struct channel_struct *dep;

	// create list of graphs for pedestals
	TList *grPedList = new TList();
	TGraphErrors *grPed;

	// create period histogram

	TString title = "Period histogram";
	TH1 *hPeriod = new TH1F(title,title, 2*((Int_t) DOMINO_NCELL), (Double_t) -DOMINO_NCELL, (Double_t) DOMINO_NCELL);

	// calculate or read pedestals from file
	grPedList = OpenPedestals(PedFile);
	grPed = (TGraphErrors *) grPedList->At(anaChannel);


	// Count number of events in data file
	int nevtDataMax = 0;
	while (!feof(fdata)) {
		fread((void *) &event_data, 1, sizeof(event_data), fdata);
		nevtDataMax++;
	}
	printf("nevtDataMax: %d\n", nevtDataMax - 1);

	if (nevt > (nevtDataMax - startEv) || nevt == 0)
		nevt = nevtDataMax - startEv;
	cout << endl << "==>> Processing " << nevt << " events from file "
			<< DataFile << endl;

	rewind(fdata);


	Int_t ievt = 1;
	// go to first event (startEv)
	while (ievt < startEv) {
		fread((void *) &event_data, 1, sizeof(event_data), fdata);
		if (feof(fdata))
			break;
		ievt++;
	}

	ievt = 1;
	Int_t flagEnd = 0;
	Int_t fitusati = 0;
	Double_t chtmp;
	Double_t PedVal, itmp;
	Double_t mean, rms;
	Double_t ratio;



	//debug canvas
	TCanvas *cfitTest = new TCanvas("cfitTest", "fit tests", 1200, 780);
	cfitTest->Divide(1,nevt);

	// loop on events

	gProgress->Reset();
	gProgress->SetMax(nevt);

	gSystem->ProcessEvents();

	while (ievt <= nevt && !flagEnd) {
		fread((void *) &event_data, 1, sizeof(event_data), fdata);
		if (feof(fdata))
			flagEnd = 1;

		p = (struct channel_struct *) &event_data.ch[0]; // read bunch of data
		dep = (struct channel_struct *) &event_data.ch[1]; // read bunch of data
		// goes to channel to analyze
		p += anaChannel;

		// read data, subtract pedestals values and save results in grAnaChDataTemp graph with
		// fixed error for each point (x = 0.5 and y = 2.1). Also generate an array with Domino
		// X and Y values

		TGraphErrors *grAnaChDataTemp = new TGraphErrors(DOMINO_NCELL);

		for (int ch = 0; ch < DOMINO_NCELL; ch++) {
			// Read pedestal value for this cell
			grPed->GetPoint(ch, itmp, PedVal);
			chtmp = (Double_t)(p->data[ch]); // data value
			chtmp = chtmp - PedVal;
			grAnaChDataTemp->SetPoint(ch, (Double_t) ch, chtmp);
			grAnaChDataTemp->SetPointError(ch, 0.5, 2.1);
		}
		// create fit functions
//.........这里部分代码省略.........