C++ TH1F::GetFunction方法代码示例

void PlotGain( Char_t *fname, Char_t *hname0, Char_t *option, UInt_t pCbc = 0, UInt_t color = 4 ){

	TH1F *h0 = new TH1F( hname0, ";V_{CTH} middle point pedestal subtracted;count", 180, -80, 0 );
	TFile *fin = new TFile( fname );
	TFile *fped = new TFile( "/Users/kirika/Work/CMS_UPGRADE/eps_calib/VCthScanCalibrationResult.root" );
	for( UInt_t i=0; i<254; i++ ){
		TString hname = Form( "h_00_%02d_%03d", pCbc, i );
		TH1F *h = (TH1F *) fin->Get( hname );
		TF1 *func = h->GetFunction( Form("f_%s", hname.Data() ) );
		h = (TH1F *) fped->Get( hname );
		TF1 *func0 = h->GetFunction( Form("f_%s", hname.Data() ) );
		h0->Fill( func->GetParameter(0) - func0->GetParameter(0) );
	}
	h0->Draw( option );
}

vector<double> one_fit(bool do_fit)
{
    TH1F* h = new TH1F("h","h", 20, -10, 10);
    h->FillRandom("gaus");
    vector<double> ret;
    if (do_fit) {
        h->Fit("gaus","");
        TF1* g = h->GetFunction("gaus");
        h->Draw();
        canvas->Modified();
        canvas->Update();
        for (int ind=0; ind < 3; ++ind) {
            ret.push_back(g->GetParameter(ind));
        }
    }
    else {
        for (int ind=0; ind < 3; ++ind) {
            ret.push_back(0.0);
        }
    }
    
    ret.push_back(h->GetMean());
    ret.push_back(h->GetRMS());
    delete h;
    return ret;
}

void PlotE( Char_t *fname, UInt_t pBeId = 0, UInt_t pFeId = 0, UInt_t pCbc = 0, UInt_t color = 4 ){


	TH1F *h0 = new TH1F("Noise", ";% ;count", 100, 0, 0.2 );
	TFile *fin = new TFile( fname );
	TFile *fped = new TFile( "/Users/kirika/Work/CMS_UPGRADE/eps_calib/VCthScanCalibrationResult.root" );
	for( UInt_t i=0; i<254; i++ ){
		TString hname = Form( "h_%02d_%02d_%02d_%03d", pBeId, pFeId, pCbc, i );
		TH1F *h = (TH1F *) fin->Get( hname );
		TF1 *func = h->GetFunction( Form("f_%s", hname.Data() ) );
		h = (TH1F *) fped->Get( hname );
		TF1 *func0 = h->GetFunction( Form("f_%s", hname.Data() ) );
		h0->Fill( func->GetParameter(1) / ( func0->GetParameter(0) - func->GetParameter(0) ) );
	}
	h0->Draw();
}

void plotMuFromToys(std::string inputFile, std::string selectString="fit_status==0"){

	// Some Global preferences
        gROOT->SetStyle("Plain");
	gSystem->Load("$CMSSW_BASE/lib/$SCRAM_ARCH/libHiggsAnalysisCombinedLimit.so");
	gStyle->SetOptFit(1111);
	gStyle->SetOptStat(1111);
	gStyle->SetPalette(1,0);

	TFile *fi_ = TFile::Open(inputFile.c_str());
        TTree *tree_sb = (TTree*) fi_->Get("tree_fit_sb");

        TH1F *mures = new TH1F("mures","",25,-2.,2.);
        
        mures->SetLineColor(kBlue+3);
        mures->SetMarkerStyle(kOpenCircle);
        mures->SetLineWidth(2);
        
        mures->GetXaxis()->SetTitle("#mu - 1");
        mures->GetYaxis()->SetTitle(Form("no toys (%d total)",int(tree_sb->GetEntries())));
        mures->GetYaxis()->SetTitleOffset(1.05);
        mures->GetXaxis()->SetTitleOffset(0.9);
        mures->GetYaxis()->SetTitleSize(0.05);
        mures->GetXaxis()->SetTitleSize(0.05);

       	TCanvas *c = new TCanvas("c","",960,800);
        tree_sb->Draw("mu>>mures",selectString.c_str());
        mures->Fit("gaus");
        mures->GetFunction("gaus")->SetLineColor(kCyan+3);
        mures->Draw("pe1");
        c->SaveAs("mlfit/mu_residual.pdf");
}

void PlotScurveGroup( Char_t *fname, UInt_t pBeId = 0, UInt_t pFeId = 0, UInt_t pCbc = 0, UInt_t pGroup = 0, UInt_t color = 4 ){

	TFile *fin = new TFile( fname );
	UInt_t cCh1(0);
	for( int i=0; i < 16; i++ ){
		cCh1 = i * 16 + pGroup*2;
		for( int j = cCh1; j < cCh1 + 2; j++ ){  
			if( j < 254 ){

				TString hname = Form( "h_%02d_%02d_%02d_%03d", pBeId, pFeId, pCbc, j );
				TH1F *h = (TH1F *) fin->Get( hname );
				h->GetXaxis()->SetRangeUser( 50, 170 );
				h->SetLineColor(1);
				h->SetMarkerColor(1);
				h->SetTitle( Form( "FE(00),CBC(%02d); VCth; Rate", pCbc ) );
				TF1 *func = h->GetFunction( Form("f_%s", hname.Data() ) );
				if( i == 0 ) h->Draw();
				else h->Draw("same");
			}
		}
	}
	TString cPadName( fname );
	cPadName.ReplaceAll( ".root", Form( "G%d.png", pGroup ) );
	gPad->Print( cPadName ); 
}

double Alignment::Fit(TH1F &beta){
  double xmin=beta.GetXaxis()->GetXmin();
  double xmax=beta.GetXaxis()->GetXmax();

  double prev_norma=-1;
  double norma,mean,sigma;


  for(int i=0;i<10;i++){
    if(xmin>beta.GetXaxis()->GetXmax() || xmax<beta.GetXaxis()->GetXmin() || xmin>=xmax) return 0;  // Fail

    //    beta.Fit("gaus","0q","",xmin,xmax);
    beta.Fit("gaus","0q","",xmin,xmax);
    TF1 *func=beta.GetFunction("gaus");
    norma=func->GetParameter(0);
    mean=func->GetParameter(1);
    sigma=fabs(func->GetParameter(2));

    xmin=mean-1.0*sigma;
    xmax=mean+2.0*sigma;
    if(fabs(norma-prev_norma)<1e-3) break;

    prev_norma=norma;

  }

  return -norma;  // This have to be minimized
}

pair <float,float> find_interpolated_point(TH2F* histo, int xbin, bool mSUGRA=true) {
  
  int minaccept=4;
  TCanvas *flatcan = new TCanvas("flatcan","flatcan");
  stringstream histoname;
  histoname << "exclusion shape for x bin " << xbin;
  TH1F *flathisto;
  if(mSUGRA) flathisto = new TH1F("flat",histoname.str().c_str(),histo->GetNbinsY(),histo->GetYaxis()->GetBinLowEdge(1),histo->GetYaxis()->GetBinLowEdge(histo->GetNbinsY())+histo->GetYaxis()->GetBinWidth(histo->GetNbinsY()));
  else flathisto = new TH1F("flat",histoname.str().c_str(),histo->GetNbinsX(),histo->GetXaxis()->GetBinLowEdge(1),histo->GetXaxis()->GetBinLowEdge(histo->GetNbinsX())+histo->GetXaxis()->GetBinWidth(histo->GetNbinsX()));
  
  int acceptedpoints=0;
  int nbins=histo->GetNbinsY();
  if(!mSUGRA) histo->GetNbinsX();
  for(int i=1;i<nbins;i++) {
    float value=0;
    if(i<=nbins-2) value=((1/3.0)*(histo->GetBinContent(xbin,i)+histo->GetBinContent(xbin+1,i)+histo->GetBinContent(xbin+2,i)));
    if(i==nbins-1) value=((1/2.0)*(histo->GetBinContent(xbin,i)+histo->GetBinContent(xbin,i+1)));
    if(i==nbins) value=(histo->GetBinContent(xbin,i));
    if(value<20&&value>0) {
      flathisto->SetBinContent(i,value);
      flathisto->SetBinError(i,TMath::Sqrt(value));
      acceptedpoints++;
    }
  }
  
  float pointone=-100;
  TLine *excluded;
  if(acceptedpoints>minaccept) {
    flathisto->Fit("expo","Q");
    TF1 *fitfunc = (TF1*)flathisto->GetFunction("expo");
    pointone=-(fitfunc->GetParameter(0)/fitfunc->GetParameter(1));
    excluded=new TLine(pointone,0,pointone,10);
  }
  
  pair <float,float> anything;
  anything.first=histo->GetXaxis()->GetBinCenter(xbin);
  anything.second=pointone;
  stringstream flatname;
  flathisto->GetYaxis()->SetRangeUser(0,10);
  flathisto->Draw();
  if(acceptedpoints>minaccept) excluded->SetLineColor(kGreen);
  if(acceptedpoints>minaccept) excluded->SetLineStyle(2);
  if(acceptedpoints>minaccept) excluded->SetLineWidth(2);
  if(acceptedpoints>minaccept) excluded->Draw("same");
  flatname << "Limits/partials/partial_" << xbin << "___" << histo->GetName() << ".png";
  if(draweachone) CompleteSave(flatcan,flatname.str());
  delete flatcan;
  delete flathisto;
  return anything;
}

//.........这里部分代码省略.........
         
         // choose start values for gen fit
         slope = (Geny.at(Geny.size()-1) - Geny.at(Geny.size()-4))/(x.at(x.size()-1) - x.at(x.size()-4));
         d_slope = slope;
         offset = Geny.at(Geny.size()-1) - (slope*x.at(x.size()-1));
         d_offset = offset;
            
         std::cout << "Gen start values: " << "slope: " << slope << " offset: " << offset << std::endl; 
         make_lin_fit(slope, d_slope, offset, d_offset);
         std::cout << "Gen fit values: " << "slope: " << slope << " offset: " << offset << std::endl; 
         
         lin_extrapol_gen->SetParameter(0, offset);
         lin_extrapol_gen->SetParError(0, d_offset);
         lin_extrapol_gen->SetParameter(1, slope);
         lin_extrapol_gen->SetParError(1, d_slope);
         extrapol_Gen->GetListOfFunctions()->Add(lin_extrapol_gen);
         
         data.reset();

         // draw extrapolations data + mc
         TCanvas *c = new TCanvas("c","",600,600);
         std::pair <float,float> minMaxPair = determineMinMax(extrapol_Data);
         c->DrawFrame(0,minMaxPair.first*0.5-0.05,alpha.back()+0.05,minMaxPair.second*1.47,(";Threshold #alpha_{max};#sigma_{A}"));
         extrapol_MC->SetMarkerStyle(20);
         extrapol_MC->SetMarkerColor(kRed+1);
         extrapol_MC->SetLineColor(kRed+1);
         extrapol_MC->Draw("P");
         extrapol_Data->SetMarkerStyle(20);
         extrapol_Data->SetMarkerColor(kBlack);
         extrapol_Data->SetLineColor(kBlack);
         extrapol_Data->Draw("Psame");
         TF1* MCTemp = new TF1();
         TF1* DataTemp = new TF1();
         extrapol_MC->GetFunction("lin_extrapol_mc")->SetLineColor(kRed+1);
         extrapol_MC->GetFunction("lin_extrapol_mc")->SetLineStyle(2);
         extrapol_Data->GetFunction("lin_extrapol_data")->SetLineColor(kBlack);
         extrapol_Data->GetFunction("lin_extrapol_data")->SetLineStyle(2);
         MCTemp=(TF1*) extrapol_MC->GetFunction("lin_extrapol_mc")->Clone();
         DataTemp=(TF1*) extrapol_Data->GetFunction("lin_extrapol_data")->Clone();
         MCTemp->SetRange(0.1,1);
         MCTemp->SetLineStyle(1);
         MCTemp->Draw("same");
         DataTemp->SetRange(0.1,1);
         DataTemp->SetLineStyle(1);
         DataTemp->Draw("same");
       
         TPaveText *label = util::LabelFactory::createPaveTextWithOffset(2,1.0,0.5);
         label->AddText("Anti-k_{T} (R=0.5) PFchs Jets");
         label->AddText( Form("%0.1f #leq |#eta| #leq %0.1f, %3.0f #leq #bar{ p}_{T} [GeV] #leq %3.0f", eta_bins[ieta], eta_bins[ieta+1], pt_bins[ipt], pt_bins[ipt+1]) );
         label->Draw("same");
  
         TLegend* leg1 = util::LabelFactory::createLegendWithOffset(2,0.6);
         leg1->AddEntry(extrapol_Data,"Extrapolation (MC smeared)","LP");
         leg1->AddEntry(extrapol_MC,"Extrapolation (MC)","LP");

         leg1->Draw();
         cmsPrel(-1, false , 8);

         TString name;
         name = Form("ClosureTest/Extrapol_Eta%i_pt%i" + suffix + ".eps", ieta, ipt);
         c->Print(name);

         TCanvas *cb = new TCanvas("c","",600,600);
         std::pair <float,float> minMaxPair2 = determineMinMax(extrapol_Data);
         cb->DrawFrame(0,minMaxPair2.first*0.5-0.05,alpha.back()+0.05,minMaxPair2.second*1.47,(";Threshold #alpha_{max, gen};#sigma_{A, gen}"));
         extrapol_Gen->SetMarkerStyle(20);

void plot1d(int hid) {
  gStyle->SetOptStat(0);
  char c[50];
  c1->Clear();  
  if(m1dHist[hid]==0){ 
    c1->Divide(2,2); 
    for(int quad=0; quad<kFgtNumQuads; quad++){
      TVirtualPad* pad1 = c1->cd(quad+1); 
      pad1->SetLogy(l1dHist[hid]);
      double xmin, xmax, ymin=0.0, ymax=0.0;
      if(l1dHist[hid]==1) ymin=0.1;
      for(int disc=0; disc<kFgtNumDiscs; disc++){
	sprintf(c,"%1d%1s-%s",disc+1,cquad[quad],c1dHist[hid]);
	//printf("Getting %s\n",c);
	TH1F *h = hist1[disc][quad][hid] = (TH1F*)file->Get(c);
	xmin=h->GetXaxis()->GetXmin();
	xmax=h->GetXaxis()->GetXmax();
	double m=h->GetMaximum();
	if(ymax<m) ymax=m;
	printf("disc=%d max=%6.1f ymax=%6.1f xmin=%6.1f xmax=%6.1f\n",disc+1,m,ymax,xmin,xmax);
      }
      sprintf(c,"Quad%1s-%s",cquad[quad],c1dHist[hid]);
      TH2F *frame = new TH2F(c,c,1,xmin,xmax,1,ymin,ymax*1.2); frame->SetStats(0); frame->Draw();
      for(int disc=0; disc<kFgtNumDiscs; disc++){
	TH1F *h=hist1[disc][quad][hid];
	h->SetLineColor(color[disc]); h->SetLineWidth(3); h->Draw("SAME");  
	if(f1dHist[hid]==0){
	  float mean=h->GetMean();
	  sprintf(c,"%1d%s mean=%6.2f",disc+1,cquad[quad],mean);
	}else if(f1dHist[hid]==1){	  
	  int res = h->Fit("gaus","0Q");
	  TF1 *f = h->GetFunction("gaus");
	  float sig = f->GetParameter(2);	
	  if(res==0 && sig>0.0001 && h->GetEntries()>5){
	    f->SetLineColor(color[disc]); f->SetLineWidth(2); f->Draw("SAME");
	    sprintf(c,"%1d%s sig=%6.3f",disc+1,cquad[quad],sig);
	  }else{
	    sprintf(c,"%1d%s",disc+1,cquad[quad]);
	  }
	}else if(f1dHist[hid]==2){	  
	  int res = h->Fit("landau","0Q");
	  TF1 *f = h->GetFunction("landau");
	  float peak = f->GetParameter(1);	
	  if(res==0 && peak>0 && h->GetEntries()>5){
	    f->SetLineColor(color[disc]); f->SetLineWidth(2); f->Draw("SAME");
	    sprintf(c,"%1d%s mpv=%6.0f",disc+1,cquad[quad],peak);
	  }else{
	    sprintf(c,"%1d%s",disc+1,cquad[quad]);
	  }
	}
	TText *t1;
	float x1= 0.2, x2= 0.55;
	float y1=0.8 - 0.07*disc;
	float y2=0.8 - 0.07*(disc-3);
        if(disc<3) { t1 = new TText(x1,y1,c); }
        else       { t1 = new TText(x2,y2,c); }
	t1->SetNDC();
	t1->SetTextSize(0.04); 
	t1->SetTextColor(color[disc]); 
	t1->Draw();
      }
    }
  }else{
    c1->Divide(4,6);
    gStyle->SetOptStat(0);
    gStyle->SetOptTitle(0);
    gStyle->SetOptFit(0);
    for(int disc=0; disc<kFgtNumDiscs; disc++){
      for(int quad=0; quad<kFgtNumQuads; quad++){
	TPad* pad = c1->cd(disc*4+quad+1);
	pad->SetRightMargin(0.01); pad->SetLeftMargin(0.1);
	pad->SetTopMargin(0.01);   pad->SetBottomMargin(0.1);
	pad->SetLogy(l1dHist[hid]);
	sprintf(c,"%1d%1s-%s",disc+1,cquad[quad],c1dHist[hid]);
	TH1F *h = hist1[disc][quad][hid] = (TH1F*)file->Get(c);
	h->SetFillColor(color[disc]);
	h->GetXaxis()->SetLabelSize(0.1);
	h->GetYaxis()->SetLabelSize(0.1);
	h->GetXaxis()->SetNdivisions(205);
	h->Draw();

	if(f1dHist[hid]==0){
	  float mean=h->GetMean();
	  sprintf(c,"%1d%s mean=%6.2f",disc+1,cquad[quad],mean);
	}else if(f1dHist[hid]==1){	  
	  int res = h->Fit("gaus","0Q");
	  TF1 *f = h->GetFunction("gaus");
	  float sig = f->GetParameter(2);	
	  if(res==0 && sig>0.0001 && h->GetEntries()>5){
	    f->SetLineColor(color[disc]); f->SetLineWidth(2); f->Draw("SAME");
	    sprintf(c,"%1d%s sig=%6.3f",disc+1,cquad[quad],sig);
	  }else{
	    sprintf(c,"%1d%s",disc+1,cquad[quad]);
	  }
	}else if(f1dHist[hid]==2){	  
	  int res = h->Fit("landau","0Q");
	  TF1 *f = h->GetFunction("landau");
	  float peak = f->GetParameter(1);	
	  if(res==0 && peak>0 && h->GetEntries()>5){
	    f->SetLineColor(1); f->SetLineWidth(2); f->Draw("SAME");
//.........这里部分代码省略.........

void fullCalib7_14(){
  const int numRuns = 6;
  const char* files[numRuns] = {
				"../crunchedFiles/run_00393_crunched.root",
				"../crunchedFiles/run_00394_crunched.root",
				"../crunchedFiles/run_00395_crunched.root",
				"../crunchedFiles/run_00398_crunched.root",
				"../crunchedFiles/run_00399_crunched.root",
				"../crunchedFiles/run_00400_crunched.root"
				};
				
  
  gStyle->SetOptFit(1);

  cout << "Which? " << endl;
  int which;
  cin >> which;


  vector<double> means;
  vector<double> sigs;
  vector<double> meanErrs;
  vector<double> sigErrs;

  for (int i = 0; i < numRuns; ++i){
    TFile* f = new TFile(files[i]);
    TTree* t = (TTree*) f->Get("t");
    cout << "got tree" << endl;
    t->Draw(Form("sipm%i.energy>>histo%i", which,i),Form("sipm%i.valid",which));
    cout << "drew" << endl;

    cout << "finding obj" << endl;
    TH1F* h = (TH1F*) gROOT->FindObject(Form("histo%i",i));
    h->Fit("gaus","0");
    means.push_back(h->GetFunction("gaus")->GetParameter(1));
    sigs.push_back(h->GetFunction("gaus")->GetParameter(2));
    meanErrs.push_back(h->GetFunction("gaus")->GetParError(1));
    sigErrs.push_back(h->GetFunction("gaus")->GetParError(2));
    cout << "end itr" << endl;
  }
  
  vector<double> variances;
  vector<double> varianceErrors;
  
  for(int i = 0; i < means.size(); ++i){
      variances.push_back(sigs[i]*sigs[i]);
      varianceErrors.push_back(2*sigs[i]*sigErrs[i]);
  }
   
    TGraphErrors* linGraph = new TGraphErrors(means.size(), &means[0], &variances[0], 
					      &meanErrs[0], &varianceErrors[0]);
    TGraphErrors* quadGraph = new TGraphErrors(means.size(), &means[0], &variances[0], 
					       &meanErrs[0], &varianceErrors[0]);
  
    linGraph->SetMarkerStyle(20);
    linGraph->Fit("pol1");
    linGraph->SetTitle(Form("SiPM %i", which));
    linGraph->GetFunction("pol1")->SetLineColor(kBlue);
    linGraph->GetXaxis()->SetTitle("Mean");
    linGraph->GetYaxis()->SetTitle("#sigma^{2}");
    quadGraph->SetMarkerStyle(20);
    quadGraph->Fit("pol2");
    quadGraph->GetFunction("pol2")->SetLineColor(kRed);
    
    linGraph->Draw("ap");
    quadGraph->Draw("psame");
  
}

void display(){
  
  gROOT->ProcessLine(".x lhcbStyle.C");
  gStyle->SetPalette(1);
	
  string filename = "../analysed/SET9THL914.root";
	TFile* file = new TFile(filename.c_str(),"OPEN");

	int nBinsX=32; int nBinsY=32;
	int nPixelsX=64; int nPixelsY=64;
	
  std::cout<<"Making global intercept tot map"<<std::endl;
  // Global intercept map (not per pixel)
	TH2F* hToTMap = new TH2F("hToTMap","hToTMap",nBinsX,0,nPixelsX,nBinsY,0,nPixelsY);
	for(int x=0;x<nBinsX;x++){
		for(int y=0;y<nBinsY;y++){
			int id = x + y*nBinsX;
			// Get the ToT distribution for this location
			string histoname = "/tpanal/ClicpixAnalysis/hClusterTOTAssociated1pixMap"+makestring(id);
			TH1F* totMapHisto = 0;
			totMapHisto = (TH1F*)gDirectory->Get(histoname.c_str());
			if(totMapHisto == 0) continue;
			if(totMapHisto->GetEntries() < 20) continue;
	
			// Fit the tot distribution and enter the value in the ToT map
			totMapHisto->Fit("gaus","0q");
			TF1* fit = totMapHisto->GetFunction("gaus");
			hToTMap->Fill(x*floor(nPixelsX/nBinsX),y*floor(nPixelsY/nBinsY),fit->GetParameter(1));
			
			delete fit; delete totMapHisto;
		}
	}
	// Plot the map of ToT
	TCanvas* canv = new TCanvas();
	hToTMap->GetXaxis()->SetTitle("Column");
	hToTMap->GetYaxis()->SetTitle("Row");
  hToTMap->GetZaxis()->SetTitle("Mean charge (ToT)");
	hToTMap->SetMaximum(9);
	hToTMap->SetMinimum(5);
	hToTMap->DrawCopy("colz");

  
  std::cout<<"Making per pixel tot map"<<std::endl;
  // Per pixel map
  TH2F* hToTMapPerPixel = new TH2F("hToTMapPerPixel","hToTMapPerPixel",nPixelsX,0,nPixelsX,nPixelsY,0,nPixelsY);
  for(int x=0;x<nPixelsX;x++){
    for(int y=0;y<nPixelsY;y++){
      int id = x + y*nPixelsX;
      // Get the ToT distribution for this location
      string histoname = "/tpanal/ClicpixAnalysis/hClusterTOTAssociated1pixPixel"+makestring(id);
      TH1F* totMapHisto = 0;
      totMapHisto = (TH1F*)gDirectory->Get(histoname.c_str());
      if(totMapHisto == 0) continue;
      if(totMapHisto->GetEntries() < 10) continue;
      
      // Fit the tot distribution and enter the value in the ToT map
      totMapHisto->Fit("gaus","0q");
      TF1* fit = totMapHisto->GetFunction("gaus");
      hToTMapPerPixel->Fill(x,y,fit->GetParameter(1));
      
      delete fit; delete totMapHisto;
    }
  }
  // Plot the map of ToT
  TCanvas* canv2 = new TCanvas();
  hToTMapPerPixel->GetXaxis()->SetTitle("Column");
  hToTMapPerPixel->GetYaxis()->SetTitle("Row");
  hToTMapPerPixel->GetZaxis()->SetTitle("Mean charge (ToT)");
  hToTMapPerPixel->SetMaximum(10);
  hToTMapPerPixel->SetMinimum(0);
  hToTMapPerPixel->DrawCopy("colz");

  std::cout<<"Making pixel response tot maps"<<std::endl;
  // Pixel response maps
  int responseBins=100;
  double responseWidth=0.05;
  TH2F* hToTMapPixelResponseX = new TH2F("hToTMapPixelResponseX","hToTMapPixelResponseX",responseBins+1,-responseWidth,responseWidth,15,0,15);
  for(double rID=0;rID<responseBins;rID++){
    // Get the ToT distribution for this location
    string histoname = "/tpanal/ClicpixAnalysis/hClusterTOTAssociated1pixResponseX"+makestringD(rID);
    TH1F* totMapHisto = 0;
    totMapHisto = (TH1F*)gDirectory->Get(histoname.c_str());
    if(totMapHisto == 0) continue;
    if(totMapHisto->GetEntries() < 1) continue;
 		// loop over all tot bins
    for(int i=0;i<45;i++){
      hToTMapPixelResponseX->Fill(2.*responseWidth*(rID/100.)-responseWidth,i,totMapHisto->GetBinContent(totMapHisto->FindBin(i)));
    }
    delete totMapHisto;
  }
  // Plot the map of ToT
  TCanvas* canv3 = new TCanvas();
  hToTMapPixelResponseX->GetXaxis()->SetTitle("Track distance x (mm)");
  hToTMapPixelResponseX->GetYaxis()->SetTitle("ToT");
  hToTMapPixelResponseX->DrawCopy("colz");

  TH1F* hToTProfilePixelResponseY = new TH1F("hToTProfilePixelResponseY","hToTProfilePixelResponseY",responseBins+1,-responseWidth,responseWidth);
  TH2F* hToTMapPixelResponseY = new TH2F("hToTMapPixelResponseY","hToTMapPixelResponseY",responseBins+1,-responseWidth,responseWidth,16,0,16);
  for(double rID=0;rID<responseBins;rID++){
    // Get the ToT distribution for this location
//.........这里部分代码省略.........

void data_to_tree()
{
	Double_t pedtime;
	Double_t pedTr;
	Double_t pedTl;
	Double_t rped;
	Double_t lped;
	//ptfactor = point to time factor
	//const Double_t ptfactor = 0.050;
	const Double_t ptfactor = 0.0407;

	TFile* f = new TFile("data_to_tree.root", "RECREATE");

	ifstream in;
	ofstream out;

	//bring data files and give name to use in this code.
	//data file name is in order type_ch0highvoltage_ch2highvoltage_#ofevent_position_3.dat
	
	const Int_t nFile = 10; // including long file for TWC 
	const Char_t* file[nFile] = {"./Cobalt_data/Co_2090v_2160v_10k_10cm_3.dat",
									"./Cobalt_data/Co_2090v_2160v_10k_20cm_3.dat",
									"./Cobalt_data/Co_2090v_2160v_10k_30cm_3.dat",
									"./Cobalt_data/Co_2090v_2160v_10k_40cm_3.dat",
									"./Cobalt_data/Co_2090v_2160v_10k_50cm_3.dat",
									"./Cobalt_data/Co_2090v_2160v_10k_60cm_3.dat",
									"./Cobalt_data/Co_2090v_2160v_10k_70cm_3.dat",
									"./Cobalt_data/Co_2090v_2160v_10k_80cm_3.dat",
									"./Cobalt_data/Co_2090v_2160v_10k_90cm_3.dat",
									"./Cobalt_data/Co_2090v_2160v_100k_50cm_3.dat"};

	const Char_t* distance[nFile] = {"10cm", "20cm", "30cm", "40cm", "50cm", "60cm", "70cm", "80cm", "90cm", "50cm100k"};

	TCanvas* c1 = new TCanvas("c1", "pedl", 200, 0, 600, 600);
	TH1F* pedlhist = new TH1F("pedlhist", "left", 200, 0, 200);
	TCanvas* c2 = new TCanvas("c2", "pedr", 210, 10, 610, 610);
	TH1F* pedrhist = new TH1F("pedrhist", "right", 200, 0, 200);

	//bring the pedestal file
	in.open("./Cobalt_data/pedestal_2090v_2160v_10k_4.dat");

	if(in.is_open())
	{
		while(1)
		{
			in >> pedtime >> pedTl >> pedTr >> lped >> rped;
			if(!in.good()) break;
			pedlhist->Fill(lped);
			pedrhist->Fill(rped);
		}
	}
	in.close();

	c1->cd();
	pedlhist->Draw();
	pedlhist->Fit("gaus", "", "", 150, 160);
	Double_t pedlmean = pedlhist->GetFunction("gaus")->GetParameter(1);
	c1->SaveAs("CFDpedl.png");

	c2->cd();
	pedrhist->Draw();
	pedrhist->Fit("gaus", "", "", 105, 125);
	Double_t pedrmean = pedrhist->GetFunction("gaus")->GetParameter(1);
	c2->SaveAs("CFDpedr.png");

	//pr: position resolution
	TTree* pr[nFile];

	for(Int_t i = 0; i < nFile; i++)
	{
		in.open(file[i]);

		Double_t time = 0;
		Double_t Tl = 0;
		Double_t Tr = 0;
		Double_t Al = 0;
		Double_t Ar = 0;

		pr[i] = new TTree(distance[i], "position resolution data");
		pr[i]->Branch("Tl", &Tl, "Tl/D");
		pr[i]->Branch("Tr", &Tr, "Tr/D");
		pr[i]->Branch("Al", &Al, "Al/D");
		pr[i]->Branch("Ar", &Ar, "Ar/D");

		//correction 
		// 1. TDCchannel * ptfactor
		// 2. ADC - pedestal mean
		if(in.is_open())
		{
			while(1)
			{
				in >> time >> Tl >> Tr >> Al >> Ar;
				if(!in.good()) break;
				Tl = Tl*ptfactor;
				Tr = Tr*ptfactor;
				Al = Al-pedlmean;
				Ar = Ar-pedrmean;
				pr[i]->Fill();
			}
		}
//.........这里部分代码省略.........

void plotTreeNorms(TTree *tree_, std::string selectString, bool do7TeV){

	// Create a map for plotting the pullsummaries:
	std::map < const char*, std::pair <double,double> > pullSummaryMap;
	int nPulls=0;

	TObjArray *l_branches = tree_->GetListOfBranches();
	int nBranches = l_branches->GetEntries();

	gStyle->SetPadTopMargin(0.01);

	TCanvas *c = new TCanvas("c","",960,800);

	std::string treename = tree_->GetName();
	c->SaveAs(Form("%s_normresiduals.pdf[",treename.c_str()));
	// File to store plots in 
	TFile *fOut = new TFile(Form("%s_normresiduals.root",treename.c_str()),"RECREATE");

        TH1F *bHd = new TH1F("bHd","",50,-1.0,1.0);
        TH1F *bHfd = new TH1F("bHfd","",50,-1.0,1.0);

	for (int iobj=0;iobj<nBranches;iobj++){

		TBranch *br =(TBranch*) l_branches->At(iobj);

		// Draw the normal histogram
		const char* name = br->GetName();

                // select only the normalizations
                string namestr(name);
                if(namestr.find("n_exp")==string::npos) continue;

                bool fitPull=true;
                bool fitPullf=true;

		double p_mean =0;

		int nToysInTree = tree_->GetEntries();
		// Find out if paramter is fitted value or constraint term.
                bool isFitted = true;
			
                p_mean = prenorms_[name].first;	// toy initial parameters from the datacards
                std::cout << "******* "<< name << " *******"<<std::endl;
                std::cout << p_mean << std::endl;
                std::cout << "******************************" <<std::endl;

                TH1F* bH = (TH1F*)bHd->Clone(Form("%s",name));
                TH1F* bHf = (TH1F*)bHfd->Clone(Form("%s_fail",name));
                
                const char* drawInput = Form("(%s-%f)/%f",name,p_mean,p_mean);
                tree_->Draw(Form("%s>>%s",drawInput,name),"");
                tree_->Draw(Form("%s>>%s_fail",drawInput,name),selectString.c_str(),"same");
                fitPull  = true;
                fitPullf = true;
                  
		bHf->SetLineColor(2);
		bH->GetXaxis()->SetTitle(bH->GetTitle());
		bH->GetYaxis()->SetTitle(Form("no toys (%d total)",nToysInTree));
		bH->GetYaxis()->SetTitleOffset(1.05);
		bH->GetXaxis()->SetTitleOffset(0.9);
		bH->GetYaxis()->SetTitleSize(0.05);
		bH->GetXaxis()->SetTitleSize(0.05);
		bH->GetXaxis()->SetTitle(Form("%s",name));
                
		
		bH->SetTitle("");	

		if ( bH->Integral() <=0 )  fitPull = false;
		if (fitPull) {bH->Fit("gaus"); bH->GetFunction("gaus")->SetLineColor(4);}
		
		if ( bHf->Integral() <=0 )  fitPullf = false;
		if (fitPullf) {bHf->Fit("gaus"); bHf->GetFunction("gaus")->SetLineColor(2);}

		c->Clear();
		
 		TPad pad1("t1","",0.01,0.01,0.66,0.95);
 		TPad pad2("t2","",0.70,0.20,0.98,0.80);

		pad1.SetNumber(1); pad2.SetNumber(2);

                if ( isFitted ) {pad2.Draw();}

		pad1.Draw();
		pad1.SetGrid(true);


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

		
		c->cd(1); bH->Draw(); bHf->Draw("same");
		TLegend *legend = new TLegend(0.6,0.8,0.9,0.89);
		legend->SetFillColor(0);
		legend->AddEntry(bH,"All Toys","L");
		legend->AddEntry(bHf,selectString.c_str(),"L");
		legend->Draw();

		if (fitPull){
			c->cd(2);
			double gap;
			TLatex *tlatex = new TLatex(); tlatex->SetNDC(); 
//.........这里部分代码省略.........

void TrigWeightFit1D(){
  TFile *f = new TFile("../data/scale_factors/HLTEffParkedABCD_1D.root");
  TFile *out = new TFile("HLTEffParkedABCD_1DFit.root","RECREATE");
  
  unsigned nVars=4;
  std::vector<std::string> histstofit[4];
  std::vector<std::string> mjjhists;
  mjjhists.push_back("hData_MJJ_1D_A");
  mjjhists.push_back("hData_MJJ_1D_BC");
  mjjhists.push_back("hData_MJJ_1D_D");
  std::vector<std::string> metl1hists;
  metl1hists.push_back("hData_MET_L1_A");
  metl1hists.push_back("hData_MET_L1_BC");
  metl1hists.push_back("hData_MET_L1_D");
  std::vector<std::string> methists;
  methists.push_back("hData_MET_1D_A");
  methists.push_back("hData_MET_1D_BC");
  methists.push_back("hData_MET_1D_D");
  std::vector<std::string> jet2hists;
  jet2hists.push_back("hData_JET2_1D_A");
  jet2hists.push_back("hData_JET2_1D_BC");
  jet2hists.push_back("hData_JET2_1D_D");

  histstofit[0]=mjjhists;
  histstofit[1]=metl1hists;
  histstofit[2]=methists;
  histstofit[3]=jet2hists;

  //create a function with 3 parameters in the range [-3,3]
  for(unsigned iVar=0;iVar<nVars;iVar++){
    for(unsigned iTurnon=0;iTurnon<histstofit[iVar].size();iTurnon++){
      TH1F *hpx = (TH1F*)f->Get(histstofit[iVar][iTurnon].c_str()); //LOOP OVER THE TRIG TURN ONS
      double rangelow,rangehigh;
      if(iVar==0){
	rangelow=600;
	rangehigh=2200;
      }
      else if(iVar==1){
	rangelow=60;
	rangehigh=400;
      }
      else if(iVar==2){
	rangelow=0;
	rangehigh=400;
      }
      else if(iVar==3){
	rangelow=0;
	rangehigh=100;
      }
      TF1 *func = new TF1("erf","0.5*(1+TMath::Erf((x-[0])/(sqrt([1]))))",rangelow,rangehigh);
      func->SetParameters(100,10000);
      func->SetParNames("Turn On","Sigma");
      func->SetParLimits(0,-10000,10000);
      func->SetParLimits(1,0,100000);
      hpx->Fit("erf","R");
      TCanvas *c1=new TCanvas();
      hpx->Draw();
      c1->SaveAs((histstofit[iVar][iTurnon]+".pdf").c_str());
      TF1 *fitresult=hpx->GetFunction("erf");
      //Get name to save as
      std::string outname=histstofit[iVar][iTurnon].substr(1);
      out->WriteObject(fitresult,("f"+outname).c_str());
    }
  }
}

void plotTree(TTree *tree_, std::string whichfit, std::string selectString){

	// Create a map for plotting the pullsummaries:
	std::map < const char*, std::pair <double,double> > pullSummaryMap;
	int nPulls=0;

	TObjArray *l_branches = tree_->GetListOfBranches();
	int nBranches = l_branches->GetEntries();

	TCanvas *c = new TCanvas("c","",960,800);

	std::string treename = tree_->GetName();
	c->SaveAs(Form("%s.pdf[",treename.c_str()));

	for (int iobj=0;iobj<nBranches;iobj++){

		TBranch *br =(TBranch*) l_branches->At(iobj);

		// Draw the normal histogram
		const char* name = br->GetName();
		bool fitPull=false;
		bool plotLH=false;

		TGraph *gr=0;
		double p_mean =0;
		double p_err  =0;

		int nToysInTree = tree_->GetEntries();
		if (doPull && findNuisancePre(name)){
			
			p_mean = bfvals_[name].first;	// toy constrainits thrown about best fit to data
			p_err  = prevals_[name].second; // uncertainties taken from card

			const char* drawInput = Form("(%s-%f)/%f",name,p_mean,p_err);
			tree_->Draw(Form("%s>>%s",drawInput,name),"");
			tree_->Draw(Form("%s>>%s_fail",drawInput,name),selectString.c_str(),"same");
			fitPull = true;
			if (doLH) {
			  gr = graphLH(name,p_err,whichfit);
			  plotLH=true;
			}
			
		}

		else{
			tree_->Draw(Form("%s>>%s",name,name),"");
			tree_->Draw(Form("%s>>%s_fail",name,name),"mu<0","same");
		}
		
		TH1F* bH  = (TH1F*) gROOT->FindObject(Form("%s",name))->Clone();
		TH1F* bHf = (TH1F*) gROOT->FindObject(Form("%s_fail",name))->Clone();
		bHf->SetLineColor(2);
		bH->GetXaxis()->SetTitle(bH->GetTitle());
		bH->GetYaxis()->SetTitle(Form("no toys (%d total)",nToysInTree));
		bH->GetYaxis()->SetTitleOffset(1.32);
		
		bH->SetTitle("");	

		if (fitPull) bH->Fit("gaus");
	
		c->Clear();
		TPad pad1("t1","",0.01,0.02,0.59,0.98);
		TPad pad2("t2","",0.59,0.04,0.98,0.62);
		TPad pad3("t3","",0.59,0.64,0.98,0.90);

		pad1.SetNumber(1); pad2.SetNumber(2); pad3.SetNumber(3);
		pad1.Draw(); pad2.Draw();pad3.Draw();
		pad2.SetGrid(true);

		c->cd(1); bH->Draw(); bHf->Draw("same");
		TLatex *titletext = new TLatex();titletext->SetNDC();titletext->SetTextSize(0.04); titletext->DrawLatex(0.1,0.95,name);
		TLegend *legend = new TLegend(0.6,0.8,0.9,0.89);
		legend->SetFillColor(0);
		legend->AddEntry(bH,"All Toys","L");
		legend->AddEntry(bHf,selectString.c_str(),"L");
		legend->Draw();

		if (doPull && plotLH) {
			c->cd(2); gr->Draw("ALP");
		}
		if (fitPull){
			c->cd(3);
			TLatex *tlatex = new TLatex(); tlatex->SetNDC(); tlatex->SetTextSize(0.12);
			tlatex->DrawLatex(0.15,0.75,Form("Mean    : %.3f #pm %.3f",bH->GetFunction("gaus")->GetParameter(1),bH->GetFunction("gaus")->GetParError(1)));
			tlatex->DrawLatex(0.15,0.60,Form("Sigma   : %.3f #pm %.3f",bH->GetFunction("gaus")->GetParameter(2),bH->GetFunction("gaus")->GetParError(2)));
			tlatex->DrawLatex(0.15,0.35,Form("Pre-fit : %.3f ",prevals_[name].first));
			tlatex->DrawLatex(0.15,0.2,Form("Best-fit (B)  : %.3f ",p_mean));
			tlatex->DrawLatex(0.15,0.05,Form("Best-fit (S+B): %.3f ",bfvals_sb_[name].first));
			
			pullSummaryMap[name]=std::make_pair<double,double>(bH->GetFunction("gaus")->GetParameter(1),bH->GetFunction("gaus")->GetParameter(2));
			nPulls++;

		}

		c->SaveAs(Form("%s.pdf",treename.c_str()));
	}
	
	if (doPull && nPulls>0){
	   
	    int nRemainingPulls = nPulls;
//.........这里部分代码省略.........