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

void DynamicExec()
{
   // Example of function called when a mouse event occurs in a pad.
   // When moving the mouse in the canvas, a second canvas shows the
   // projection along X of the bin corresponding to the Y position
   // of the mouse. The resulting histogram is fitted with a gaussian.
   // A "dynamic" line shows the current bin position in Y.
   // This more elaborated example can be used as a starting point
   // to develop more powerful interactive applications exploiting CINT
   // as a development engine.
   //
   // Author:  Rene Brun
   
   TObject *select = gPad->GetSelected();
   if(!select) return;
   if (!select->InheritsFrom("TH2")) {gPad->SetUniqueID(0); return;}
   TH2 *h = (TH2*)select;
   gPad->GetCanvas()->FeedbackMode(kTRUE);

   //erase old position and draw a line at current position
   int pyold = gPad->GetUniqueID();
   int px = gPad->GetEventX();
   int py = gPad->GetEventY();
   float uxmin = gPad->GetUxmin();
   float uxmax = gPad->GetUxmax();
   int pxmin = gPad->XtoAbsPixel(uxmin);
   int pxmax = gPad->XtoAbsPixel(uxmax);
   if(pyold) gVirtualX->DrawLine(pxmin,pyold,pxmax,pyold);
   gVirtualX->DrawLine(pxmin,py,pxmax,py);
   gPad->SetUniqueID(py);
   Float_t upy = gPad->AbsPixeltoY(py);
   Float_t y = gPad->PadtoY(upy);

   //create or set the new canvas c2
   TVirtualPad *padsav = gPad;
   TCanvas *c2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c2");
   if(c2) delete c2->GetPrimitive("Projection");
   else   c2 = new TCanvas("c2","Projection Canvas",710,10,700,500);
   c2->SetGrid();
   c2->cd();

   //draw slice corresponding to mouse position
   Int_t biny = h->GetYaxis()->FindBin(y);
   TH1D *hp = h->ProjectionX("",biny,biny);
   hp->SetFillColor(38);
   char title[80];
   sprintf(title,"Projection of biny=%d",biny);
   hp->SetName("Projection");
   hp->SetTitle(title);
   hp->Fit("gaus","ql");
   hp->GetFunction("gaus")->SetLineColor(kRed);
   hp->GetFunction("gaus")->SetLineWidth(6);
   c2->Update();
   padsav->cd();
}

void GetCosmicOut(){
  gSystem->Load("lib/libtest.so");
  IDHandler* handler  = new IDHandler("Data/crystal.txt");
  CsIImage*  image    = new CsIImage(handler);
  image->SetTitle("Distribution On image; x[mm]; y[mm]");
  TH1D* gainHist      = new TH1D("gain","Distribution of gain;Integrated ADC;N/10cnt",500,0,5000);				 
  
  std::ofstream ofs("CH_Cosmic_outList_3599_3600.txt");  
  
  Double_t gain;
  Double_t Sigma;
  Double_t Norm;
  Double_t HistEntries;
  TFile* tf = new TFile("cosmic_3599_3600.root");
  
  for( int i = 0; i< 2716; i++){    
    TH1D* hist  = (TH1D*)tf->Get(Form("his_CH%04d",i));
    TF1*  func  = NULL;
    func        = hist->GetFunction("landau");
    if( func == NULL ){
      continue;
    }
    if( i < 2240){      
      gain = func->GetParameter(1);      
      image->Fill(i, func->GetParameter(1));
      ofs << i << "\t" 
	  << gain << "\t"
	  << func->GetParameter(2) << std::endl; 
    }else{
      gain = func->GetParameter(1)/2;
      image->Fill(i, gain);
      ofs << i << "\t" 
	  << gain << "\t" 
	  << func->GetParameter(2)/2   << std::endl; 
    }
    gainHist->Fill(gain);
  }

  ofs.close();
  TCanvas* can = new TCanvas("can","",1600,800);
  can->Divide(2,1);  
  can->cd(1);
  image->Draw();
  can->cd(2);
  gainHist->Draw();

}

// Show the Pt spectra, and start the tree viewer.
int showEventSample()
{
  // Load needed libraries
  loadLibraries();

  // Open the file
  TFile* file = TFile::Open(FILENAME, "READ");
  if (!file || !file->IsOpen()) {
    Error("showEventSample", "Couldn;t open file %s", FILENAME);
    return 1;
  }

  // Get the tree
  TTree* tree = (TTree*)file->Get(TREENAME);
  if (!tree) {
    Error("showEventSample", "couldn't get TTree %s", TREENAME);
    return 2;
  }

  // Start the viewer.
  tree->StartViewer();

  // Get the histogram
  TH1D* hist = (TH1D*)file->Get(HISTNAME);
  if (!hist) {
    Error("showEventSample", "couldn't get TH1D %s", HISTNAME);
    return 4;
  }

  // Draw the histogram in a canvas
  gStyle->SetOptStat(1);
  TCanvas* canvas = new TCanvas("canvas", "canvas");
  canvas->SetLogy();
  hist->Draw("e1");
  TF1* func = hist->GetFunction("expo");

  char expression[64];
  sprintf(expression,"T #approx %5.1f", -1000 / func->GetParameter(1));
  TLatex* latex = new TLatex(1.5, 1e-4, expression);
  latex->SetTextSize(.1);
  latex->SetTextColor(4);
  latex->Draw();

  return 0;
}

//.........这里部分代码省略.........
  tp4->Draw();
  tp5->Draw();
  tp6->Draw();
  tp7->Draw();
  tp8->Draw();
  tp9->Draw();
  tp10->Draw();
  tp11->Draw();
  tp12->Draw();
  tp13->Draw();
  tp14->Draw();
  tp15->Draw();
  tp16->Draw();
  tp17->Draw();
  tp18->Draw();

  ic = 2;
  Int_t ipad=0;
  for(Int_t i=1; i<4; i++){
    TString hisname = "h_xy_module_L"; hisname  += i;
    hisname += "t_halfmodule_axial_sensor0";
    TString layer = "L"; layer+=i; layer+= "TA"; 
    TH2D *his2 = (TH2D*) f->Get(hisname.Data());
    TH1D *his = (TH1D*)his2->ProjectionY();
    ipad++;
    TString pd = "tp"; pd+=ipad;
    TPad *pad = (TPad*) cTop->GetListOfPrimitives()->FindObject(pd.Data());
    pad->cd();
    Double_t low = his->GetMean()-w*his->GetRMS();
    Double_t up = his->GetMean()+w*his->GetRMS();    
    his->Fit("gaus","Q0","",low,up);
    his2->Draw("colz");
    module->Draw("same");
    TF1 *fit = his->GetFunction("gaus");
    // cout << "sensor " << i << " axial mean :  " << his->GetMean()*1000. << " - RMS : " << his->GetRMS()*1000. << " (um) " << " " << his->GetEntries() << endl;
    // cout << "sensor " << i << " axial mu :  " << fit->GetParameter(1)*1000. << " - sigma : " << fit->GetParameter(2)*1000. << " (um) " << endl;
    // outf << layer.Data() << " " << fit->GetParameter(1)*1000. << " " << fit->GetParameter(2)*1000. << " " << his->GetEntries() << endl;

    hisname = "h_xy_module_L"; hisname  += i;
    hisname += "t_halfmodule_stereo_sensor0";
    layer = "L"; layer+=i; layer+= "TS"; 
    his2 = (TH2D*) f->Get(hisname.Data());
    his = (TH1D*)his2->ProjectionY();
    ipad++;
    pd = "tp"; pd+=ipad;
    pad = (TPad*) cTop->GetListOfPrimitives()->FindObject(pd.Data());
    pad->cd();
    low = his->GetMean()-w*his->GetRMS();
    up = his->GetMean()+w*his->GetRMS();    
    his->Fit("gaus","Q","",low,up); his2->Draw("colz");
    module->Draw("same");
    fit = his->GetFunction("gaus");
  }
  //  ic = 19;
  for(Int_t i=4; i<7; i++){
    TString hisname = "h_xy_module_L"; hisname  += i;
    TString hisname2 = hisname;
    hisname += "t_halfmodule_axial_hole_sensor0";
    hisname2 += "t_halfmodule_axial_slot_sensor0";
    TH2D *hisOther = (TH2D*) f->Get(hisname2.Data());
    Double_t norma2 = hisOther->GetMaximum();
    TString layer = "L"; layer+=i; layer+= "TAHo"; 
    TH2D *his2 = (TH2D*) f->Get(hisname.Data());
    TH1D *his = (TH1D*)his2->ProjectionY();     
    Double_t norma1 = his2->GetMaximum();
    Double_t normamax = TMath::Max(norma1,norma2);

//.........这里部分代码省略.........
  leg->SetTextFont(42);
  leg->SetFillStyle(0);
  leg->AddEntry(h,"Data","pl");
  leg->AddEntry(f,"Fit","l");
  leg->AddEntry(mass,"D^{0}+#bar{D^{#lower[0.2]{0}}} Signal","f");
  leg->AddEntry(massSwap,"K-#pi swapped","f");
  leg->AddEntry(background,"Combinatorial","l");
  leg->Draw("same");

  TLatex Tl;
  Tl.SetNDC();
  Tl.SetTextAlign(12);
  Tl.SetTextSize(0.04);
  Tl.SetTextFont(42);
  Tl.DrawLatex(0.18,0.93, "#scale[1.25]{CMS} Performance");
  Tl.DrawLatex(0.65,0.93, Form("%s #sqrt{s_{NN}} = 5.02 TeV",collisionsystem.Data()));

  TLatex* tex;

  tex = new TLatex(0.22,0.78,Form("%.1f < p_{T} < %.1f GeV/c",ptmin,ptmax));
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.04);
  tex->SetLineWidth(2);
  tex->Draw();

  tex = new TLatex(0.22,0.83,"|y| < 1.0");
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.04);
  tex->SetLineWidth(2);
  tex->Draw();

  h->GetFunction(Form("f%d",count))->Delete();
  TH1F* histo_copy_nofitfun = ( TH1F * ) h->Clone("histo_copy_nofitfun");
  histo_copy_nofitfun->Draw("esame");

 TH1D* hTest = new TH1D("hTest","",nbinsmasshisto,minhisto,maxhisto);

 for (int m=0;m<yieldtotal;m++){
 double r = f->GetRandom();
 hTest->Fill(r);
 
 }

 TF1* ffaketotal=(TF1*)f->Clone("ffake");
 TF1* ffakemass=(TF1*)mass->Clone("ffakemass");
 TF1* ffakebackground=(TF1*)background->Clone("ffakebackground");
 TF1* ffakemassSwap=(TF1*)massSwap->Clone("ffakemassSwap");
 
 Double_t yieldtotal_original = ffaketotal->Integral(minhisto,maxhisto)/binwidthmass;
 Double_t yieldmass_original = ffakemass->Integral(minhisto,maxhisto)/binwidthmass;
 Double_t yieldbackground_original = ffakebackground->Integral(minhisto,maxhisto)/binwidthmass;
 Double_t yieldswapped_original = ffakemassSwap->Integral(minhisto,maxhisto)/binwidthmass;

 TH1D* hTestFake = new TH1D("hTestFake","",nbinsmasshisto,minhisto,maxhisto);
  ffakemass->SetParameter(2,ffaketotal->GetParameter(2)*0.8);
  ffakemass->SetParameter(10,ffaketotal->GetParameter(10)*0.8);
 
Double_t yieldmass_modified= ffakemass->Integral(minhisto,maxhisto)/binwidthmass;

cout<<"mass original="<<yieldmass_original<<endl;
cout<<"mass modified="<<yieldmass_modified<<endl;

 for (int m=0;m<yieldmass_original*scalefactor;m++){
   double r = ffakemass->GetRandom();

// nEvents is how many events we want.
int makeEventSample(Int_t nEvents)
{
  // Load needed libraries
  loadLibraries();

  // Create an instance of the Pythia event generator ...
  TPythia6* pythia = new TPythia6;

  // ... and initialise it to run p+p at sqrt(200) GeV in CMS
  pythia->Initialize("cms", "p", "p", 200);

  // Open an output file
  TFile* file = TFile::Open(FILENAME, "RECREATE");
  if (!file || !file->IsOpen()) {
    Error("makeEventSample", "Couldn;t open file %s", FILENAME);
    return 1;
  }

  // Make a tree in that file ...
  TTree* tree = new TTree(TREENAME, "Pythia 6 tree");

  // ... and register a the cache of pythia on a branch (It's a
  // TClonesArray of TMCParticle objects. )
  TClonesArray* particles = (TClonesArray*)pythia->GetListOfParticles();
  tree->Branch(BRANCHNAME, &particles);

  // Now we make some events
  for (Int_t i = 0; i < nEvents; i++) {
    // Show how far we got every 100'th event.
    if (i % 100 == 0)
      cout << "Event # " << i << endl;

    // Make one event.
    pythia->GenerateEvent();

    // Maybe you want to have another branch with global event
    // information.  In that case, you should process that here.
    // You can also filter out particles here if you want.

    // Now we're ready to fill the tree, and the event is over.
    tree->Fill();
  }

  // Show tree structure
  tree->Print();

  // After the run is over, we may want to do some summary plots:
  TH1D* hist = new TH1D(HISTNAME, "p_{#perp}  spectrum for  #pi^{+}",
                        100, 0, 3);
  hist->SetXTitle("p_{#perp}");
  hist->SetYTitle("dN/dp_{#perp}");
  char expression[64];
  sprintf(expression,"sqrt(pow(%s.fPx,2)+pow(%s.fPy,2))>>%s",
          BRANCHNAME, BRANCHNAME, HISTNAME);
  char selection[64];
  sprintf(selection,"%s.fKF==%d", BRANCHNAME, PDGNUMBER);
  tree->Draw(expression,selection);

  // Normalise to the number of events, and the bin sizes.
  hist->Sumw2();
  hist->Scale(3 / 100. / hist->Integral());
  hist->Fit("expo", "QO+", "", .25, 1.75);
  TF1* func = hist->GetFunction("expo");
  func->SetParNames("A", "- 1 / T");
  // and now we flush and close the file
  file->Write();
  file->Close();

  return 0;
}

//.........这里部分代码省略.........
    if (hphac!=0) {hphac->SetFillColor(0);hphac->SetLineColor(4);hphac->SetMarkerColor(4);hphac->Draw("same");lgd3->AddEntry(hphac,"p_{had}");}
    lgd3->Draw();
    
    TCanvas *c2 = new TCanvas();
    TLegend *lgd2 = new TLegend(0.15, 0.6,0.4,0.85);
    if(h2pka!=0){h2pka->SetFillColor(0);h2pka->SetLineColor(1);h2pka->SetMarkerColor(1);h2pka->Draw();lgd2->AddEntry(h2pka,"p_{K}");} else continue;
    if(h2pee!=0){h2pee->SetFillColor(0);h2pee->SetLineColor(2);h2pee->SetMarkerColor(2);h2pee->Draw("same");lgd2->AddEntry(h2pee,"p_{e}");}
    if(h2pmu!=0){h2pmu->SetFillColor(0);h2pmu->SetLineColor(3);h2pmu->SetMarkerColor(3);h2pmu->Draw("same");lgd2->AddEntry(h2pmu,"p_{#mu}");}
    if(h2pha!=0){h2pha->SetFillColor(0);h2pha->SetLineColor(4);h2pha->SetMarkerColor(4);h2pha->Draw("same");lgd2->AddEntry(h2pha,"p_{had}");}
    lgd2->Draw();
    
    TCanvas *c4 = new TCanvas();
    TLegend *lgd4 = new TLegend(0.15, 0.6,0.4,0.85);
    if (h2pkac!=0) {h2pkac->SetFillColor(0);h2pkac->SetLineColor(1);h2pkac->SetMarkerColor(1);h2pkac->Draw(); lgd4->AddEntry(h2pkac,"p_{K}");} else continue;
    if (h2peec!=0) {h2peec->SetFillColor(0);h2peec->SetLineColor(2);h2peec->SetMarkerColor(2);h2peec->Draw("same");lgd4->AddEntry(h2peec,"p_{e}");}
    if (h2pmuc!=0) {h2pmuc->SetFillColor(0);h2pmuc->SetLineColor(3);h2pmuc->SetMarkerColor(3);h2pmuc->Draw("same");lgd4->AddEntry(h2pmuc,"p_{#mu}");}
    if (h2phac!=0) {h2phac->SetFillColor(0);h2phac->SetLineColor(4);h2phac->SetMarkerColor(4);h2phac->Draw("same");lgd4->AddEntry(h2phac,"p_{had}");}
    lgd4->Draw();
    
    sprintf(filename,"%s_1p",file[enei]);
    outfile->WriteTObject(c1,filename);
    sprintf(filename,"%s_2p",file[enei]);
    outfile->WriteTObject(c2,filename);
    sprintf(filename,"%s_1p_cut",file[enei]);
    outfile->WriteTObject(c3,filename);
    sprintf(filename,"%s_2p_cut",file[enei]);
    outfile->WriteTObject(c4,filename);

    double p_ka = sqrt(pow(ene[enei]/2,2)-mka*mka);
    double p_mu = sqrt(pow(ene[enei]/2,2)-mmu*mmu);
    //std::cout<<p_ka<<' '<<p_mu<<std::endl;
    hpka->Fit("gaus","R","",p_ka-0.1,p_ka+0.1);
    hpmu->Fit("gaus","R","",p_mu-0.1,p_mu+0.1);
    pka[enei] = hpka->GetFunction("gaus")->GetParameter(1);
    pkae[enei]= hpka->GetFunction("gaus")->GetParameter(2);
    pka3e[enei]= 3*pkae[enei]/2.6*3;
    pka5e[enei]= 5*pkae[enei]/2.6*3;
    pmu[enei] = hpmu->GetFunction("gaus")->GetParameter(1);
    pmue[enei]= hpmu->GetFunction("gaus")->GetParameter(2);
    pmu3e[enei]= 3*pmue[enei]/2.6*3;
    pmu5e[enei]= 5*pmue[enei]/2.6*3;
    std::cout<<p_ka<<' ' << pka[enei]<<" "<< pkae[enei]<<std::endl;
    std::cout<<p_mu<<' ' << pmu[enei]<<" "<< pmue[enei]<<std::endl;
    
    delete filei1;
    delete filei2;
    delete filei3;
    delete filei4;
    delete c1;
    delete c2;
    delete c3;
    delete c4;
    delete legend;
    delete lgd2;
    delete lgd3;
    delete lgd4;
  }

  new TCanvas();
  TGraphErrors *graph1 = new TGraphErrors(21,ene,pka,enee,pkae);
  graph1->SetFillColor(0);
  graph1->SetMarkerStyle(24);
  TGraphErrors *graph1_1 = new TGraphErrors(21,ene,pka,enee,pkae);
  graph1_1->SetFillColor(3);
  TGraphErrors *graph1_3 = new TGraphErrors(21,ene,pka,enee,pka3e);
  graph1_3->SetFillColor(4);

// ----------------------------------------------------------------------
void invertedMuonID(string era = "2016BF", float bdt = -99., TH1D *hresult = 0) {
  if ("all" == era) {
    invertedMuonID("2016BF");
    invertedMuonID("2016GH");
    return;
  }

  TFile *fd(0), *fm0(0), *fm1(0);
  float bdtcut(0.);
  cout << "======================================================================" << endl;
  if (era == "2016BF") {
    fd = TFile::Open("/scratch/ursl/bmm4/s01/bmm-data-bmmCharmonium2016BF-s01.root");
    fm0 = TFile::Open("/scratch/ursl/bmm4/s01/bmm-mc-RunIISpring16DR80-BdToMuMu-2016BF-s01.root");
    fm1 = TFile::Open("/scratch/ursl/bmm4/s01/bmm-mc-RunIISpring16DR80-BsToMuMu-2016BF-s01.root");
    if (bdt < -1.) {
      bdtcut = 0.30;
    } else {
      bdtcut = bdt;
    }
  } else if (era == "2016GH") {
    fd = TFile::Open("/scratch/ursl/bmm4/s01/bmm-data-bmmCharmonium2016GH-s01.root");
    fm0 = TFile::Open("/scratch/ursl/bmm4/s01/bmm-mc-RunIISpring16DR80-BdToMuMu-2016GH-s01.root");
    fm1 = TFile::Open("/scratch/ursl/bmm4/s01/bmm-mc-RunIISpring16DR80-BsToMuMu-2016GH-s01.root");
    if (bdt < -1.) {
      bdtcut = 0.31;
    } else {
      bdtcut = bdt;
    }
  } else if (era == "2012") {
    fd = TFile::Open("/scratch/ursl/bmm4/s01/bmm-data-bmmMuOnia2012-s01.root");
    fm0 = TFile::Open("/scratch/ursl/bmm4/s01/bmm-mc-Winter17_private-BdToMuMu-s01.root");
    fm1 = TFile::Open("/scratch/ursl/bmm4/s01/bmm-mc-Winter17_private-BsToMuMu-s01.root");
    if (bdt < -1.) {
      bdtcut = 0.34;
    } else {
      bdtcut = bdt;
    }
  } else if (era == "2011") {
    fd = TFile::Open("/scratch/ursl/bmm4/s01/bmm-data-bmmMuOnia2011-s01.root");
    fm0 = TFile::Open("/scratch/ursl/bmm4/s01/bmm-mc-Summer17_private-BdToMuMu-s01.root");
    fm1 = TFile::Open("/scratch/ursl/bmm4/s01/bmm-mc-Summer17_private-BsToMuMu-s01.root");
    if (bdt < -1.) {
      bdtcut = 0.28;
    } else {
      bdtcut = bdt;
    }
  }

  TTree *td = (TTree*)fd->Get("candAnaMuMu/events");
  TTree *tm0 = (TTree*)fm0->Get("candAnaMuMu/events");
  TTree *tm1 = (TTree*)fm1->Get("candAnaMuMu/events");

  TH1D *hd0 = new TH1D("hd0", "mass", 60, 4.8, 6.0);
  TH1D *hm0 = new TH1D("hm0", "mass", 60, 4.8, 6.0);
  TH1D *hm1 = new TH1D("hm1", "mass", 60, 4.8, 6.0);
  string cut(Form("chan==0 && !gmuid &&hlt1&&tos&&l1t &&bdt>%4.3f && m < 6 && m > 4.8", bdtcut));
  cout << cut << endl;
  td->Draw("m>>hd0", cut.c_str());
  tm0->Draw("m>>hm0", cut.c_str());
  tm1->Draw("m>>hm1", cut.c_str());

  c0->Clear();
  // c0->Divide(2,2);
  // c0->cd(1);
  hd0->Draw("hist");
  // B0 -> mu mu
  // c0->cd(2);
  hm0->Fit("gaus");
  double p0 = hm0->GetFunction("gaus")->GetParameter(1);
  double s0 = hm0->GetFunction("gaus")->GetParameter(2);
  cout << "p0: " << p0 << " s0: " << s0 << endl;
  // Bs -> mu mu
  // c0->cd(3);
  hm1->Fit("gaus");
  double p1 = hm1->GetFunction("gaus")->GetParameter(1);
  double s1 = hm1->GetFunction("gaus")->GetParameter(2);
  cout << "p1: " << p1 << " s1: " << s1 << endl;

  // c0->cd(4);
  TH1D *hFit = (TH1D*)hd0->Clone("hFit");
  hFit->SetTitle("");
  hFit->SetMinimum(0.);
  TF1 *f1 = new TF1("f1", adhocfunction, 4.8, 6.0, 7);
  f1->SetParameter(0, 1);  f1->SetParLimits(0, 0., 1000.);
  f1->SetParameter(1, p0); f1->FixParameter(1, p0);
  f1->SetParameter(2, s0); f1->FixParameter(2, s0);
  f1->SetParameter(3, 1);  f1->SetParLimits(3, 0., 1000.);
  f1->SetParameter(4, p1); f1->FixParameter(4, p1);
  f1->SetParameter(5, s1); f1->FixParameter(5, s1);
  f1->SetParameter(6, 1);

  hFit->Fit(f1, "RL", "", 5.1, 6.0);

  double nb0 = hFit->GetFunction("f1")->GetParameter(0);
  double eb0 = hFit->GetFunction("f1")->GetParError(0);
  double nbs = hFit->GetFunction("f1")->GetParameter(3);
  double ebs = hFit->GetFunction("f1")->GetParError(3);

  TF1 *f0 = (TF1*)hm0->GetFunction("gaus")->Clone("f0");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  mass->SetParError(4,f->GetParError(9));
  mass->SetParError(5,f->GetParError(10));
  mass->SetFillColor(kOrange-3);
  mass->SetFillStyle(3002);
  mass->SetLineColor(kOrange-3);
  mass->SetLineWidth(3);
  mass->SetLineStyle(2);
  
  TF1* massSwap = new TF1(Form("fmassSwap%d",count),"[0]*(1-[2])*Gaus(x,[1],[3])/(sqrt(2*3.14159)*[3])");
  massSwap->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(7),f->GetParameter(8));
  massSwap->SetParError(0,f->GetParError(0));
  massSwap->SetParError(1,f->GetParError(1));
  massSwap->SetParError(2,f->GetParError(7));
  massSwap->SetParError(3,f->GetParError(8));
  massSwap->SetFillColor(kGreen+4);
  massSwap->SetFillStyle(3005);
  massSwap->SetLineColor(kGreen+4);
  massSwap->SetLineWidth(3);
  massSwap->SetLineStyle(1);
  
  h->SetXTitle("m_{#piK} (GeV/c^{2})");
  h->SetYTitle("Entries / (5 MeV/c^{2})");
  h->GetXaxis()->CenterTitle();
  h->GetYaxis()->CenterTitle();
  h->SetAxisRange(0,h->GetMaximum()*1.4*1.2,"Y");
  h->GetXaxis()->SetTitleOffset(1.3);
  h->GetYaxis()->SetTitleOffset(1.8);
  h->GetXaxis()->SetLabelOffset(0.007);
  h->GetYaxis()->SetLabelOffset(0.007);
  h->GetXaxis()->SetTitleSize(0.045);
  h->GetYaxis()->SetTitleSize(0.045);
  h->GetXaxis()->SetTitleFont(42);
  h->GetYaxis()->SetTitleFont(42);
  h->GetXaxis()->SetLabelFont(42);
  h->GetYaxis()->SetLabelFont(42);
  h->GetXaxis()->SetLabelSize(0.04);
  h->GetYaxis()->SetLabelSize(0.04);
  h->SetMarkerSize(0.8);
  h->SetMarkerStyle(20);
  h->SetStats(0);
  h->Draw("e");

  background->Draw("same");   
  mass->SetRange(minhisto,maxhisto);	
  mass->Draw("same");
  massSwap->SetRange(minhisto,maxhisto);
  massSwap->Draw("same");
  f->Draw("same");
  
  Double_t yield = mass->Integral(minhisto,maxhisto)/binwidthmass;
  Double_t yieldErr = mass->Integral(minhisto,maxhisto)/binwidthmass*mass->GetParError(0)/mass->GetParameter(0);
  
  std::cout<<"YIELD="<<yield<<"and error"<<yieldErr<<std::endl;
  std::cout<<"relative error="<<yieldErr/yield<<std::endl;

  TLegend* leg = new TLegend(0.65,0.58,0.82,0.88,NULL,"brNDC");
  leg->SetBorderSize(0);
  leg->SetTextSize(0.04);
  leg->SetTextFont(42);
  leg->SetFillStyle(0);
  leg->AddEntry(h,"Data","pl");
  leg->AddEntry(f,"Fit","l");
  leg->AddEntry(mass,"D^{0}+#bar{D^{#lower[0.2]{0}}} Signal","f");
  leg->AddEntry(massSwap,"K-#pi swapped","f");
  leg->AddEntry(background,"Combinatorial","l");
  leg->Draw("same");

  TLatex Tl;
  Tl.SetNDC();
  Tl.SetTextAlign(12);
  Tl.SetTextSize(0.04);
  Tl.SetTextFont(42);
  Tl.DrawLatex(0.18,0.93, "#scale[1.25]{CMS} Preliminary");
  Tl.DrawLatex(0.65,0.93, Form("%s #sqrt{s_{NN}} = 5.02 TeV",collisionsystem.Data()));

  TLatex* tex;

  tex = new TLatex(0.22,0.78,Form("%.1f < %s  < %.1f",ptmin,variableplot.Data(),ptmax));
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.04);
  tex->SetLineWidth(2);
  tex->Draw();

  tex = new TLatex(0.22,0.83,"|y| < 1.0");
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.04);
  tex->SetLineWidth(2);
  tex->Draw();

  h->GetFunction(Form("f%d",count))->Delete();
  TH1F* histo_copy_nofitfun = ( TH1F * ) h->Clone("histo_copy_nofitfun");
  histo_copy_nofitfun->Draw("esame");
//
  if(nBins==1) c->SaveAs(Form("DMass-inclusive%s_%d.pdf",collisionsystem.Data(),count));
  else c->SaveAs(Form("DMass%s_%d.pdf",collisionsystem.Data(),count));
  
  return mass;
}

//.........这里部分代码省略.........
		}
		
		hAut->SetBinError(ibin, binError);
		
	}//end Asym calc
	
	
	
	//DRAW HISTOGRAMS
	
	hInvarM->Draw();
	
	TCanvas* cNup = new TCanvas();
	hNumberUp->Draw();
	
	TCanvas* cNdown = new TCanvas();
	hNumberDown->Draw();
	
	TCanvas* cAut = new TCanvas();
	cAut->SetName("cAut");
	
	TF1* fitFunc = new TF1("fitFunc","[0]*cos(x)+[1]",-PI,PI);
	hAut->Fit("fitFunc","R");
	hAut->Draw();
	
	hAut->GetXaxis()->SetTitle("#phi_{S} - #phi_{R}");
	
	char title[150];
	sprintf(title, "%.1f < P_{T}^{#pi^{1}#pi^{2}} < %.1f  %.1f < M_{inv}^{#pi^{1}#pi^{2}} < %.1f  %.1f < #eta^{#pi^{1}#pi^{2}} < %.1f", ptCutLo, ptCutHi, mCutLo, mCutHi, etaCutLo, etaCutHi);
	
	hAut->SetTitle(title);
	
	
	cout << "Asym  = " << hAut->GetFunction("fitFunc")->GetParameter(0) << endl;
	cout << "error = " << hAut->GetFunction("fitFunc")->GetParError(0) << endl;
	cout << "chi2  = " << hAut->GetFunction("fitFunc")->GetChisquare() << endl;
	cout << "Nup   = " << hNumberUp->GetEntries() << endl;
	cout << "Ndown = " << hNumberDown->GetEntries() << endl;
  cout << "ran5  = " << rand5 << endl;
  cout << "ran6  = " << rand6 << endl;
  cout << "ran9  = " << rand9 << endl;
  cout << "ran10 = " << rand10 << endl;
	cout << "blueU = " << blueU << endl;
	cout << "blueD = " << blueD << endl;
	cout << "yellowU = " << yellowU << endl;
	cout << "yellowD = " << yellowD << endl;
	
	
	
	
	//SAVE CANVAS	
	
	
	stringstream pt;
	stringstream eta;
	stringstream mass;
	
	string part1 = "_ptBin";
	string part2 = "_massBin";
	string part3 = "_etaBin";
	
	string ptBinStr;
	string etaBinStr;
	string massBinStr;
	
	pt << ptBin;

void readPullBias(char*prefix="pull_bias"){

  gSystem->mkdir(Form("figures_%s",prefix));
  char* dec[2] = {"EB","EE"};
  const int fBinsPt[]={21,23,26,30,35,40,45,50,60,85,120,300};
  const int nPtBin = sizeof(fBinsPt)/sizeof(fBinsPt[0])-1;
  const int NRUNS=8;

  ofstream foutm;
  foutm.open("mean.dat");

  ofstream fouts;
  fouts.open("sigma.dat");

  TH1D* hpull = new TH1D("hpull","",50,-5.0,5.0);
  hpull->SetXTitle("(Fit-Input)/Fit_err");
  TH1D* hbias = new TH1D("hbias","",100,-0.5,0.5);
  hbias->SetXTitle("(Fit-Input)/Input");

  TCanvas* c1 = new TCanvas("c1","",500,500);

  char tmp[1000];
  for(int ieta=0; ieta < 2; ieta++){
    for(int ipt=0; ipt < nPtBin; ipt++){
      hpull->Reset();
      hbias->Reset();
      for(int irun=0; irun < NRUNS; irun++){
	sprintf(tmp,"%s%d_%s_pt%d.dat",prefix,irun,dec[ieta],fBinsPt[ipt]);		cout << "Opening file " << tmp << endl;
	FILE *infile =  fopen(tmp,"r");  
	cout << "Opened file " << tmp << endl;

	Double_t pull=-100., bias=-100.; // combined isolation, pt, eta
	int flag = 1;
	if(infile==NULL){
	  cout << "File does not exist" << endl;
	  continue;
	}
	  
	while (flag!=-1){
	  flag =fscanf(infile,"%lf %lf",&pull, &bias);
	  if(flag!=-1){
	    hpull->Fill(pull);
	    hbias->Fill(bias);
	  }
 
	}// keep reading files as long as text exists

      } // end of iterations
      
      cout << hpull->GetEntries() << endl;
      cout << hbias->GetEntries() << endl;
 
      hpull->SetTitle(Form("%s, pt = %d-%d GeV",dec[ieta],
			   fBinsPt[ipt],fBinsPt[ipt+1]));
      
      hbias->SetTitle(Form("%s, pt = %d-%d GeV",dec[ieta],
			   fBinsPt[ipt],fBinsPt[ipt+1]));
      if(hpull->GetEntries()<1)continue;
      hpull->Fit("gaus");
      TF1* f1 = hpull->GetFunction("gaus");
      foutm << dec[ieta] << " pt" << fBinsPt[ipt] << " pull: " << f1->GetParameter(1) << " +- " << f1->GetParError(1) << endl;
      fouts << dec[ieta] << " pt" << fBinsPt[ipt] << " pull: " << f1->GetParameter(2) << " +- " << f1->GetParError(2) << endl;

      c1->Print(Form("figures_%s/pull_%s_%s_pt%d.png",
		     prefix,prefix,dec[ieta],fBinsPt[ipt]));

      if(hbias->GetEntries()<1)continue;
      hbias->Fit("gaus");
      TF1* f2 = hbias->GetFunction("gaus");
      foutm << dec[ieta] << " pt" << fBinsPt[ipt] << " bias: " << f2->GetParameter(1) << " +- " << f2->GetParError(1) << endl;
      fouts << dec[ieta] << " pt" << fBinsPt[ipt] << " bias: " << f2->GetParameter(2) << " +- " << f2->GetParError(2) << endl;
      c1->Print(Form("figures_%s/bias_%s_%s_pt%d.png",
		     prefix,prefix,dec[ieta],fBinsPt[ipt]));

    } // end of loop over pt bins
  } // end of loop over eta bins



}

//.........这里部分代码省略.........
  TF1 *gaussPlusLine = new TF1("gaussPlusLine", "[0] + [1]*x + [2]*TMath::Gaus(x,[3],[4],1)");

  bkgMEUS->Fit(line, "", "", 2.1, 2.5);
  double offset = line->GetParameter(0);
  double slope = line->GetParameter(1);
  const double LcMass = 2.28646;

  gaussPlusLine->FixParameter(0,offset);
  gaussPlusLine->FixParameter(1,slope);
  gaussPlusLine->SetParameter(3,LcMass);

  gaussPlusLine->SetParameter(2,1);
  gaussPlusLine->SetParameter(4,sigma);

  sig->Fit(gaussPlusLine, "", "", 2.1, 2.5);

  // -- Drawing
  sig->GetXaxis()->SetTitle("#font[12]{m}_{pK#pi} (GeV/#font[12]{c}^{2})");
  sig->GetXaxis()->CenterTitle();
  sig->GetXaxis()->SetLabelSize(0.04);
  sig->GetXaxis()->SetTitleSize(0.055);
  sig->GetXaxis()->SetTitleOffset(0.8);
  // sig->GetXaxis()->SetTickLength(0.001);
  sig->GetYaxis()->SetTitle("#font[12]{N}/(10 MeV/#font[12]{c}^{2})");
  sig->GetYaxis()->CenterTitle();
  sig->GetYaxis()->SetLabelSize(0.04);
  sig->GetYaxis()->SetTitleSize(0.055);
  sig->GetYaxis()->SetTitleOffset(0.8);
  // sig->GetYaxis()->SetRangeUser(17, 95);

  sig->SetMarkerStyle(20);
  sig->SetLineColor(kBlack);

  bkgMEUS->SetMarkerStyle(21);
  bkgMEUS->SetMarkerColor(kOrange+5);
  bkgMEUS->SetLineColor(kOrange+5);
  bkgMEUS->GetFunction("line")->SetLineColor(kRed);
  bkgMEUS->GetFunction("line")->SetLineStyle(2);
  bkgMEUS->GetFunction("line")->SetNpx(10000);

  bkgSELS->SetMarkerStyle(22);
  bkgSELS->SetMarkerColor(kBlue-7);
  bkgSELS->SetLineColor(kBlue-7);

  bkgMELS->SetMarkerStyle(23);
  bkgMELS->SetMarkerColor(kGreen+2);
  bkgMELS->SetLineColor(kGreen+2);
  
  sig->GetFunction("gaussPlusLine")->SetLineColor(kRed);
  sig->GetFunction("gaussPlusLine")->SetNpx(10000);

  sig->GetXaxis()->SetRangeUser(2.1,2.5);
  sig->Draw("E");
  bkgMEUS->Draw("Esame");
  bkgSELS->Draw("Esame");
  // bkgMELS->Draw("Esame");

  TLegend *leg = new TLegend(0.5,0.6,0.89,0.8);
  leg->SetFillStyle(0);
  leg->SetBorderSize(0);
  leg->AddEntry(sig, "Same event correct sign", "pl");
  leg->AddEntry(bkgSELS, "Same event wrong sign", "pl");
  leg->AddEntry(bkgMEUS, "Mixed event correct sign", "pl");
  // leg->AddEntry(bkgMELS, "Mixed event wrong sign", "pl");
  leg->Draw();
  // sideBandMarker->Draw("same BAR");

  TLegend *dataSet = new TLegend(0., 0.7, 0.4, 0.89);
  dataSet->SetFillStyle(0);
  dataSet->SetBorderSize(0);
  dataSet->AddEntry("", "#font[22]{Au+Au 200 GeV, 10-60%}","");
  dataSet->AddEntry("",Form("%1.0f GeV/#font[12]{c} < #font[12]{p}_{T} < %1.0f GeV/#font[12]{c}", minPt, maxPt),"");
  dataSet->Draw();

  // --  yield calculation
  const float min = 2.26;
  const float max = 2.31;

  // -- get part of the Gauss curve
  TF1 *normGaus = new TF1("normGaus", "TMath::Gaus(x,[0],[1],1)");
  double mean = gaussPlusLine->GetParameter(3);
  double width= gaussPlusLine->GetParameter(4);
  normGaus->SetParameter(0,mean);
  normGaus->SetParameter(1,width);
  double norm = normGaus->Integral(min, max);

  int minBin = sig->FindBin(min) + 1;
  int maxBin = sig->FindBin(max);
  cout << "integral from " << sig->GetBinCenter(minBin) - 0.005 << " to " << sig->GetBinCenter(maxBin) + 0.005 << endl;
  float NSig = sig->Integral(minBin,maxBin);
  float NBkg = bkgMEUS->Integral(minBin, maxBin);
  float Nlc = (NSig - NBkg)/norm;
  float error = std::sqrt(NSig/norm);

  cout << "N Sig = " << NSig << endl;
  cout << "N Bkg = " << NBkg << endl;
  cout << "norm = " << norm << endl;
  cout << "N Lambda_c = " << Nlc << endl;
  cout << "sigma = " << error << endl;
}

void ABCD2vari_p3(const char* ISO) { //Only iso is external, all the regions go in a loop!!! to be fed in the macro at once

  
  TString PREFIX = "./";

  //Get the fractions
  //gStyle->SetOptStat(0);
  TFile* fB = new TFile(PREFIX+"factor_qcd_"+TString(ISO)+"_B.root", "open");
  fB->cd();
  TH1D* hRatioB = (TH1D*)gDirectory->Get("hfactor");
  hRatioB->Fit("pol0");
  //2D needs, preliminary
  TH1D* hRatioB2030 = (TH1D*)gDirectory->Get("hfactor2030");
  hRatioB2030->Fit("pol0");
  TH1D* hRatioB3045 = (TH1D*)gDirectory->Get("hfactor3045");
  hRatioB3045->Fit("pol0");
  TH1D* hRatioB4560 = (TH1D*)gDirectory->Get("hfactor4560");
  hRatioB4560->Fit("pol0");
  TH1D* hRatioB60120 = (TH1D*)gDirectory->Get("hfactor60120");
  hRatioB60120->Fit("pol0");
  TH1D* hRatioB120200 = (TH1D*)gDirectory->Get("hfactor120200");
  hRatioB120200->Fit("pol0");
  TH1D* hRatioB2001500 = (TH1D*)gDirectory->Get("hfactor2001500");
  hRatioB2001500->Fit("pol0");

  TF1 *myfitB = (TF1*)hRatioB->GetFunction("pol0");
  Double_t factor_B = myfitB->GetParameter(0);
  Double_t factor_B_err = myfitB->GetParError(0);
  //2D needs, preliminary
  TF1 *myfitB2030 = (TF1*)hRatioB2030->GetFunction("pol0");
  Double_t factor_B2030 = myfitB2030->GetParameter(0);
  Double_t factor_B_err2030 = myfitB2030->GetParError(0);
  TF1 *myfitB3045 = (TF1*)hRatioB3045->GetFunction("pol0");
  Double_t factor_B3045 = myfitB3045->GetParameter(0);
  Double_t factor_B_err3045 = myfitB3045->GetParError(0);
  TF1 *myfitB4560 = (TF1*)hRatioB4560->GetFunction("pol0");
  Double_t factor_B4560 = myfitB4560->GetParameter(0);
  Double_t factor_B_err4560 = myfitB4560->GetParError(0);
  TF1 *myfitB60120 = (TF1*)hRatioB60120->GetFunction("pol0");
  Double_t factor_B60120 = myfitB60120->GetParameter(0);
  Double_t factor_B_err60120 = myfitB60120->GetParError(0);
  TF1 *myfitB120200 = (TF1*)hRatioB120200->GetFunction("pol0");
  Double_t factor_B120200 = myfitB120200->GetParameter(0);
  Double_t factor_B_err120200 = myfitB120200->GetParError(0);
  TF1 *myfitB2001500 = (TF1*)hRatioB2001500->GetFunction("pol0");
  Double_t factor_B2001500 = myfitB2001500->GetParameter(0);
  Double_t factor_B_err2001500 = myfitB2001500->GetParError(0);


  //repeat for region D
  TFile* fD = new TFile(PREFIX+"factor_qcd_"+TString(ISO)+"_D.root", "open");
  fD->cd();
  TH1D* hRatioD = (TH1D*)gDirectory->Get("hfactor");
  hRatioD->Fit("pol0");

  //2D needs, preliminary
  TH1D* hRatioD2030 = (TH1D*)gDirectory->Get("hfactor2030");
  hRatioD2030->Fit("pol0");
  TH1D* hRatioD3045 = (TH1D*)gDirectory->Get("hfactor3045");
  hRatioD3045->Fit("pol0");
  TH1D* hRatioD4560 = (TH1D*)gDirectory->Get("hfactor4560");
  hRatioD4560->Fit("pol0");
  TH1D* hRatioD60120 = (TH1D*)gDirectory->Get("hfactor60120");
  hRatioD60120->Fit("pol0");
  TH1D* hRatioD120200 = (TH1D*)gDirectory->Get("hfactor120200");
  hRatioD120200->Fit("pol0");
  TH1D* hRatioD2001500 = (TH1D*)gDirectory->Get("hfactor2001500");
  hRatioD2001500->Fit("pol0");

  TF1 *myfitD = (TF1*)hRatioD->GetFunction("pol0");
  Double_t factor_D = myfitD->GetParameter(0);
  Double_t factor_D_err = myfitD->GetParError(0);
  //2D needs, preliminary
  TF1 *myfitD2030 = (TF1*)hRatioD2030->GetFunction("pol0");
  Double_t factor_D2030 = myfitD2030->GetParameter(0);
  Double_t factor_D_err2030 = myfitD2030->GetParError(0);
  TF1 *myfitD3045 = (TF1*)hRatioD3045->GetFunction("pol0");
  Double_t factor_D3045 = myfitD3045->GetParameter(0);
  Double_t factor_D_err3045 = myfitB3045->GetParError(0);
  TF1 *myfitD4560 = (TF1*)hRatioD4560->GetFunction("pol0");
  Double_t factor_D4560 = myfitD4560->GetParameter(0);
  Double_t factor_D_err4560 = myfitD4560->GetParError(0);
  TF1 *myfitD60120 = (TF1*)hRatioD60120->GetFunction("pol0");
  Double_t factor_D60120 = myfitD60120->GetParameter(0);
  Double_t factor_D_err60120 = myfitD60120->GetParError(0);
  TF1 *myfitD120200 = (TF1*)hRatioD120200->GetFunction("pol0");
  Double_t factor_D120200 = myfitD120200->GetParameter(0);
  Double_t factor_D_err120200 = myfitD120200->GetParError(0);
  TF1 *myfitD2001500 = (TF1*)hRatioD2001500->GetFunction("pol0");
  Double_t factor_D2001500 = myfitD2001500->GetParameter(0);
  Double_t factor_D_err2001500 = myfitD2001500->GetParError(0);

  // directory with data
  TString protocol = "file://";
  TString dirname = "/scratch/lustreC/a/asvyatko/DY2013/rootfiles/";

  // EWK
  TFileCollection* c6 = new TFileCollection("WJets","WJets");
  c6->Add(protocol+dirname+"WJets_PU"+"/*.root");

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

//.........这里部分代码省略.........
		fMix->SetParLimits(0,0,2*integralTotalYield);
		fMix->SetParLimits(1,0,2*integralTotalYield);

		fMix->SetLineColor(2);
		fMix->SetFillColor(kRed-9);
		fMix->SetFillStyle(1001);

		float fitRangeL = 0;
		float fitRangeH = 0.08;

		hD0DcaData->GetXaxis()->SetRangeUser(0,0.07);
		hD0DcaData->Draw();
		int fitStatus = 1;
		TFitResultPtr fitResult;
		double fitPrecision = 1.e-6;
		while(fitStatus)
		{
			TFitter::SetPrecision(fitPrecision);
			fMix->SetParameters(0.5*integralTotalYield,0.5*integralTotalYield);
			fMix->SetParError(0,0.1*integralTotalYield);
			fMix->SetParError(1,0.1*integralTotalYield);
			fitResult = hD0DcaData->Fit("fMix","E SNQ0", "", fitRangeL, fitRangeH);
			fitStatus = fitResult->Status();
			cout<<"fit precision: "<<TFitter::GetPrecision()<<"   status: "<<fitStatus<<endl;
			if(fitStatus)
				fitPrecision *= 10;
		}
		cout<<"============== do main fit ============"<<endl;
		fMix->SetParameters(integralTotalYield,0.9);
		fMix->SetParError(0,0.1*integralTotalYield);
		fMix->SetParError(1,0.1);
		fMix->SetNpx(10000);
		fitResult = hD0DcaData->Fit("fMix","E S0", "", fitRangeL, fitRangeH);
		hD0DcaData->GetFunction("fMix")->Draw("flsame");
		fitStatus = fitResult->Status();
		cout<<"fit precision: "<<TFitter::GetPrecision()<<"   status: "<<fitStatus<<endl;

		TF1* fNP = new TF1("fNP",&funNonPrompt, 0., 0.5, 2);
		fNP->SetParameters(fMix->GetParameter(0),fMix->GetParameter(1));
		fNP->SetRange(fitRangeL,fitRangeH);
		fNP->SetLineColor(4);
		fNP->SetFillStyle(1001);
		fNP->SetFillColor(kBlue-9);
		fNP->SetNpx(10000);
		fNP->Draw("same");  

		hD0DcaData->Draw("same");

		promptDYield[i-1] = fMix->GetParameter(0);
		promptDYieldErrorDataOnly[i-1] = fMix->GetParError(0);
		bToDYield[i-1] = fMix->GetParameter(1);
		bToDYieldErrorDataOnly[i-1] = fMix->GetParError(1);
		totalYield[i-1] = promptDYield[i-1]+bToDYield[i-1];
		promptFraction[i-1] = promptDYield[i-1]/totalYield[i-1];

		cout<<"chi2 / NDF: "<<fitResult->Chi2()<<" / "<<fitResult->Ndf()<<endl;

		texCms->Draw();
		texCol->Draw();
		texPtY->Draw();

		TLatex* texPrompt = new TLatex(0.4,0.73,Form("Prompt D^{0} yield : %.0f #pm %.0f",fMix->GetParameter(0),fMix->GetParError(0)));
		texPrompt->SetNDC();
		texPrompt->SetTextFont(42);
		texPrompt->SetTextSize(0.06);
		texPrompt->SetLineWidth(2);

//.........这里部分代码省略.........
	      if (wx[5][iS]!=0) {
		if (iL>22) x = 10*(2*wx[4][iS]+wx[3][iS]-wx[1][iS]-2*wx[0][iS])/wx[5][iS];
		else x = 10*(wx[2][iS]-wx[0][iS])/wx[5][iS];
	      }
	      double y = 0;//10*(wy[2][iS]-wy[0][iS])/wy[3][iS];
	      if (wy[5][iS]!=0) {
		if (iL>22) y = 10*(2*wy[4][iS]+wy[3][iS]-wy[1][iS]-2*wy[0][iS])/wy[5][iS];
		else y = 10*(wy[2][iS]-wy[0][iS])/wy[5][iS];
	      }

	      //if (fabs(y-yt)>5) std::cout << " --- iL=" << iL << " iS=" << iS 
	      //<< " x=" << x << " xt=" << xt 
	      //<< " y=" << y << " yt=" << yt 
	      //<< std::endl;
	      p_posx[ipu][iL][iS]->Fill(x-xt);
	      p_posy[ipu][iL][iS]->Fill(y-yt);
	    }

	  }//loop on layers
	  if (savePoint) p_intercalibSigmaSquare[ipu]->Fill(Etotsq);
	}//loop on entries
    
      }//loop on pt

      //fill first point with linear weighting
      TLatex lat;
      char buf[500];

      for (unsigned iL(0);iL<nLayers;++iL){
	mycFit->cd();
	TH1D *projy = p_deltavsreco_x[ipu][iL]->ProjectionY();
	projy->Draw();
	projy->Fit("gaus","0+Q");
	TF1 *fitx = projy->GetFunction("gaus");
	fitx->SetLineColor(6);
	fitx->Draw("same");
	sprintf(buf,"Layer %d, linear weighting, pu=%d",iL,pu[ipu]);
	lat.DrawLatexNDC(0.1,0.96,buf);
	
	grX[ipu][iL]->SetPoint(0,0.5,fitx->GetParameter(2));
	grX[ipu][iL]->SetPointError(0,0,fitx->GetParError(2));
	grXrms[ipu][iL]->SetPoint(0,0.5,projy->GetRMS());
	grXrms[ipu][iL]->SetPointError(0,0,projy->GetRMSError());
	mycFit->Update();
	if (savePoint) mycFit->Print("PLOTS/fits_x.pdf");
	
	mycFit->cd();
	projy = p_deltavsreco_y[ipu][iL]->ProjectionY();
	projy->Draw();
	projy->Fit("gaus","0+Q");
	TF1 *fity = projy->GetFunction("gaus");
	fitx->SetLineColor(6);
	fitx->Draw("same");
	sprintf(buf,"Layer %d, linear weighting, pu=%d",iL,pu[ipu]);
	lat.DrawLatexNDC(0.1,0.96,buf);
	
	grY[ipu][iL]->SetPoint(0,0.5,fity->GetParameter(2));
	grY[ipu][iL]->SetPointError(0,0,fity->GetParError(2));
	grYrms[ipu][iL]->SetPoint(0,0.5,projy->GetRMS());
	grYrms[ipu][iL]->SetPointError(0,0,projy->GetRMSError());
	mycFit->Update();
	if (savePoint) mycFit->Print("PLOTS/fits_y.pdf");
      }

      //fit vs w0, get w0min
      for (unsigned iL(0);iL<nLayers;++iL){