C++ TF1::DrawCopy方法代码示例

void mathBeta() {
  TCanvas *c1=new TCanvas("c1", "TMath::BetaDist",600,800);
  c1->Divide(1, 2);
  TVirtualPad *pad1 = c1->cd(1);
  pad1->SetGrid();
  TF1 *fbeta = new TF1("fbeta", "TMath::BetaDist(x, [0], [1])", 0, 1);
  fbeta->SetParameters(0.5, 0.5);
  TF1 *f1 = fbeta->DrawCopy();
  f1->SetLineColor(kRed);
  f1->SetLineWidth(1);
  fbeta->SetParameters(0.5, 2);
  TF1 *f2 = fbeta->DrawCopy("same");
  f2->SetLineColor(kGreen);
  f2->SetLineWidth(1);
  fbeta->SetParameters(2, 0.5);
  TF1 *f3 = fbeta->DrawCopy("same");
  f3->SetLineColor(kBlue);
  f3->SetLineWidth(1);
  fbeta->SetParameters(2, 2);
  TF1 *f4 = fbeta->DrawCopy("same");
  f4->SetLineColor(kMagenta);
  f4->SetLineWidth(1);
  TLegend *legend1 = new TLegend(.5,.7,.8,.9);
  legend1->AddEntry(f1,"p=0.5  q=0.5","l");
  legend1->AddEntry(f2,"p=0.5  q=2","l");
  legend1->AddEntry(f3,"p=2    q=0.5","l");
  legend1->AddEntry(f4,"p=2    q=2","l");
  legend1->Draw();

  TVirtualPad *pad2 = c1->cd(2);
  pad2->SetGrid();
  TF1 *fbetai=new TF1("fbetai", "TMath::BetaDistI(x, [0], [1])", 0, 1);
  fbetai->SetParameters(0.5, 0.5);
  TF1 *g1=fbetai->DrawCopy();
  g1->SetLineColor(kRed);
  g1->SetLineWidth(1);
  fbetai->SetParameters(0.5, 2);
  TF1 *g2=fbetai->DrawCopy("same");
  g2->SetLineColor(kGreen);
  g2->SetLineWidth(1);
  fbetai->SetParameters(2, 0.5);
  TF1 *g3=fbetai->DrawCopy("same");
  g3->SetLineColor(kBlue);
  g3->SetLineWidth(1);
  fbetai->SetParameters(2, 2);
  TF1 *g4=fbetai->DrawCopy("same");
  g4->SetLineColor(kMagenta);
  g4->SetLineWidth(1);

  TLegend *legend2 = new TLegend(.7,.15,0.9,.35);
  legend2->AddEntry(f1,"p=0.5  q=0.5","l");
  legend2->AddEntry(f2,"p=0.5  q=2","l");
  legend2->AddEntry(f3,"p=2    q=0.5","l");
  legend2->AddEntry(f4,"p=2    q=2","l");
  legend2->Draw();
  c1->cd();
}

void mathLaplace(){
   TCanvas *c1=new TCanvas("c1", "TMath::LaplaceDist",600,800);
   c1->Divide(1, 2);
   TVirtualPad *pad1 = c1->cd(1);
   pad1->SetGrid();
   TF1 *flaplace = new TF1("flaplace", "TMath::LaplaceDist(x, [0], [1])", -10, 10);
   flaplace->SetParameters(0, 1);
   TF1 *f1 = flaplace->DrawCopy();
   f1->SetLineColor(kRed);
   f1->SetLineWidth(1);
   flaplace->SetParameters(0, 2);
   TF1 *f2 = flaplace->DrawCopy("same");
   f2->SetLineColor(kGreen);
   f2->SetLineWidth(1);
   flaplace->SetParameters(2, 1);
   TF1 *f3 = flaplace->DrawCopy("same");
   f3->SetLineColor(kBlue);
   f3->SetLineWidth(1);
   flaplace->SetParameters(2, 2);
   TF1 *f4 = flaplace->DrawCopy("same");
   f4->SetLineColor(kMagenta);
   f4->SetLineWidth(1);
   TLegend *legend1 = new TLegend(.7,.7,.9,.9);
   legend1->AddEntry(f1,"alpha=0 beta=1","l");
   legend1->AddEntry(f2,"alpha=0 beta=2","l");
   legend1->AddEntry(f3,"alpha=2 beta=1","l");
   legend1->AddEntry(f4,"alpha=2 beta=2","l");
   legend1->Draw();

   TVirtualPad *pad2 = c1->cd(2);
   pad2->SetGrid();
   TF1 *flaplacei=new TF1("flaplacei", "TMath::LaplaceDistI(x, [0], [1])", -10, 10);
   flaplacei->SetParameters(0, 1);
   TF1 *g1=flaplacei->DrawCopy();
   g1->SetLineColor(kRed);
   g1->SetLineWidth(1);
   flaplacei->SetParameters(0, 2);
   TF1 *g2=flaplacei->DrawCopy("same");
   g2->SetLineColor(kGreen);
   g2->SetLineWidth(1);
   flaplacei->SetParameters(2, 1);
   TF1 *g3=flaplacei->DrawCopy("same");
   g3->SetLineColor(kBlue);
   g3->SetLineWidth(1);
   flaplacei->SetParameters(2, 2);
   TF1 *g4=flaplacei->DrawCopy("same");
   g4->SetLineColor(kMagenta);
   g4->SetLineWidth(1);

   TLegend *legend2 = new TLegend(.7,.15,0.9,.35);
   legend2->AddEntry(f1,"alpha=0 beta=1","l");
   legend2->AddEntry(f2,"alpha=0 beta=2","l");
   legend2->AddEntry(f3,"alpha=2 beta=1","l");
   legend2->AddEntry(f4,"alpha=2 beta=2","l");
   legend2->Draw();
   c1->cd();
}

// tutorial illustrating the use of TMath::GammaDist and TMath::LogNormal
void mathGammaNormal(){
  TCanvas *myc = new TCanvas("c1","gamma and lognormal",10,10,600,800);
  myc->Divide(1,2);
  TPad *pad1 = myc->cd(1);
  pad1->SetLogy();
  pad1->SetGrid();

  //TMath::GammaDist
  TF1 *fgamma = new TF1("fgamma", "TMath::GammaDist(x, [0], [1], [2])", 0, 10);
  fgamma->SetParameters(0.5, 0, 1);
  TF1 *f1 = fgamma->DrawCopy();
  f1->SetMinimum(1e-5);
  f1->SetLineColor(kRed);
  fgamma->SetParameters(1, 0, 1);
  TF1 *f2 = fgamma->DrawCopy("same");
  f2->SetLineColor(kGreen);
  fgamma->SetParameters(2, 0, 1);
  TF1 *f3 = fgamma->DrawCopy("same");
  f3->SetLineColor(kBlue);
  fgamma->SetParameters(5, 0, 1);
  TF1 *f4 = fgamma->DrawCopy("same");
  f4->SetLineColor(kMagenta);
  TLegend *legend1 = new TLegend(.2,.15,.5,.4);
  legend1->AddEntry(f1,"gamma = 0.5 mu = 0  beta = 1","l");
  legend1->AddEntry(f2,"gamma = 1   mu = 0  beta = 1","l");
  legend1->AddEntry(f3,"gamma = 2   mu = 0  beta = 1","l");
  legend1->AddEntry(f4,"gamma = 5   mu = 0  beta = 1","l");
  legend1->Draw();

  //TMath::LogNormal
  TPad *pad2 = myc->cd(2);
  pad2->SetLogy();
  pad2->SetGrid();
  TF1 *flog = new TF1("flog", "TMath::LogNormal(x, [0], [1], [2])", 0, 5);
  flog->SetParameters(0.5, 0, 1);
  TF1 *g1 = flog->DrawCopy();
  g1->SetLineColor(kRed);
  flog->SetParameters(1, 0, 1);
  TF1 *g2 = flog->DrawCopy("same");
  g2->SetLineColor(kGreen);
  flog->SetParameters(2, 0, 1);
  TF1 *g3 = flog->DrawCopy("same");
  g3->SetLineColor(kBlue);
  flog->SetParameters(5, 0, 1);
  TF1 *g4 = flog->DrawCopy("same");
  g4->SetLineColor(kMagenta);
  TLegend *legend2 = new TLegend(.2,.15,.5,.4);
  legend2->AddEntry(g1,"sigma = 0.5 theta = 0  m = 1","l");
  legend2->AddEntry(g2,"sigma = 1   theta = 0  m = 1","l");
  legend2->AddEntry(g3,"sigma = 2   theta = 0  m = 1","l");
  legend2->AddEntry(g4,"sigma = 5   theta = 0  m = 1","l");
  legend2->Draw();
}

void tStudent()
{

   //gSystem->Load("libMathMore");
   // this is the way to force load of MathMore in Cling
   ROOT::Math::MathMoreLibrary::Load();

   int n=100;
   double a=-5.;
   double b=5.;
   //double r  = 3;
   TF1* pdf = new TF1("pdf", "ROOT::Math::tdistribution_pdf(x,3.0)", a,b);
   TF1* cum = new TF1("cum", "ROOT::Math::tdistribution_cdf(x,3.0)", a,b);

   TH1D* quant = new TH1D("quant", "", 9, 0, 0.9);

   for(int i=1; i < 10; i++)
      quant->Fill((i-0.5)/10.0, ROOT::Math::tdistribution_quantile((1.0*i)/10, 3.0 ) );

   double xx[10];
   xx[0] = -1.5;
   for(int i=1; i<9; i++)
      xx[i]= quant->GetBinContent(i);
   xx[9] = 1.5;
   TH1D* pdfq[10];
   //int nbin = n/10.0;
   for(int i=0; i < 9; i++) {
      int nbin = n * (xx[i+1]-xx[i])/3.0 + 1.0;
      TString name = "pdf";
      name += i;
      pdfq[i]= new TH1D(name, "", nbin,xx[i],xx[i+1] );
      for(int j=1; j<nbin; j++) {
         double x= j*(xx[i+1]-xx[i])/nbin + xx[i];
         pdfq[i]->SetBinContent(j, ROOT::Math::tdistribution_pdf(x,3));
      }
   }

   TCanvas *Canvas = new TCanvas("DistCanvas", "Student Distribution graphs", 10, 10, 800, 700);
   pdf->SetTitle("Student t distribution function");
   cum->SetTitle("Cumulative for Student t");
   quant->SetTitle("10-quantiles  for Student t");
   Canvas->Divide(2, 2);
   Canvas->cd(1);
   pdf->SetLineWidth(2);
   pdf->DrawCopy();
   Canvas->cd(2);
   cum->SetLineWidth(2);
   cum->SetLineColor(kRed);
   cum->Draw();
   Canvas->cd(3);
   quant->Draw();
   quant->SetLineWidth(2);
   quant->SetLineColor(kBlue);
   quant->SetStats(0);
   Canvas->cd(4);
   pdfq[0]->SetTitle("Student t & its quantiles");
   pdf->SetTitle("");
   pdf->Draw();
   //pdfq[0]->SetAxisRange(-1.5, 0, 1.5,1.0);
   pdfq[0]->SetTitle("Student t & its quantiles");
   for(int i=0; i < 9; i++) {
      pdfq[i]->SetStats(0);
      pdfq[i]->SetFillColor(i+1);
      pdfq[i]->Draw("same");
   }
   Canvas->Modified();
   Canvas->cd();
}

//.........这里部分代码省略.........
    if (doNorm) {
      systHists120[bdt][1]->Scale(GetXsection(115.0)*GetBR(115.0));
      systHists120[bdt][2]->Scale(GetXsection(125.0)*GetBR(125.0));
      systHists135[bdt][1]->Scale(GetXsection(130.0)*GetBR(130.0));
      systHists135[bdt][2]->Scale(GetXsection(140.0)*GetBR(140.0));
    }
    if (doNorm) {
      int120->Scale(norm120/int120->Integral());
      int135->Scale(norm135/int135->Integral());
    }
    
    TList *int120list = new TList();
    int120list->Add(systHists120[bdt][1]);
    int120list->Add(systHists120[bdt][2]);
    int120list->Add(int120);
    plotFrac(int120list,systHists120[bdt][0],"syst120",false);
    
    TList *int135list = new TList();
    int135list->Add(systHists135[bdt][1]);
    int135list->Add(systHists135[bdt][2]);
    int135list->Add(int135);
    plotFrac(int135list,systHists135[bdt][0],"syst135",false);
    TH1F *frac=(TH1F*)systHists135[bdt][0]->Clone();
    frac->Add(int135,-1);
    frac->Divide(int135);
    TH1F *linFrac = (TH1F*)linearBin(frac);
    linFrac->GetYaxis()->SetRangeUser(-0.4,0.4);
    TCanvas *c = new TCanvas();
    systTestF->cd();
    linFrac->Write();
    linFrac->Draw();
    TF1 *fit = new TF1("fit","[0]*(x-36.)^2",0.,36.);
    linFrac->Fit(fit,"q");
    fit->DrawCopy("same");
    fit->SetLineColor(kGray);
    fit->SetLineStyle(2);
    double var = fit->GetParameter(0);
    double err = fit->GetParError(0);
    fit->SetParameter(0,var+err);
    fit->DrawCopy("same");
    fit->SetParameter(0,var-err);
    fit->DrawCopy("same");

    c->Print(Form("plots/%s/fracs/systTest135.png",BDTtype[bdt].c_str()),"png");

    // get plus and minus templates
    std::pair<TH1F*,TH1F*> result135 = GetPMSyst(systHists135[bdt][0],int135,"th1f_sig_"+BDTtype[bdt]+"_ggh_135.0_cat0_sigInt");
    TH1F* sigIntDown = (TH1F*)linearBin(result135.first);
    TH1F* sigIntUp = (TH1F*)linearBin(result135.second);
    TH1F* sigIntCent = (TH1F*)linearBin(int135);
    diagFile << (result135.first)->GetName() << std::endl;
    diagFile << (result135.second)->GetName() << std::endl;
    
    TCanvas *test = new TCanvas();
    sigIntDown->SetLineColor(2);
    sigIntUp->SetLineColor(4);
    sigIntDown->Draw();
    sigIntUp->Draw("same");
    sigIntCent->Draw("same'");
    test->Print(("plots/"+BDTtype[bdt]+"/fracs/PMsysts135.png").c_str(),"png");
   
    outFile->cd();
    (result135.first)->Write();
    (result135.second)->Write();
   
  }

//.........这里部分代码省略.........
  TH1D* ppSysR = new TH1D("ppSysR","",1,0,1);
  ppSysR->SetBinContent(1,ppSysR60);

  
  
  // xjg distributions
  TCanvas *c1 = new TCanvas("c1","",1100,330);
  makeMultiPanelCanvas(c1,4,1,0.0,0.0,0.24,0.15,0.075);
  c1->cd(0);
  drawCMSppPbPb(0.1,0.95);
  //c1->Divide(4,1,0.0,0.0);
  for ( int icent=0; icent<=3 ; icent++) {
    c1->cd( 4 - icent);
    //hxgj[khimc][icent]->SetAxisRange(-.2,2.5,"Y");
    hxgj[khimc][icent]->SetAxisRange(0,2.5,"Y");
    hxgj[khimc][icent]->SetNdivisions(505);
    //    hxgj[khimc][icent]->SetTitle(";x_{J#gamma} = p^{Jet}_{T}/p^{#gamma}_{T}; #frac{1}{N_{J#gamma}} #frac{dN_{J#gamma}}{dx_{J#gamma}}");
    hxgj[khimc][icent]->SetTitle(";x_{J#gamma}; #frac{1}{N_{J#gamma}} #frac{dN_{J#gamma}}{dx_{J#gamma}}");
    handsomeTH1(hxgj[khimc][icent]);
    fixedFontHist(hxgj[khimc][icent],1,1.35);
    mcStyle2(hxgj[khimc][icent]);
    handsomeTH1(hxgj[khidata][icent],2);
    
    hxgj[khimc][icent]->GetYaxis()->SetTitleOffset(1.5);
    TH1D * htemp41 = (TH1D*)hxgj[khimc][icent]->Clone(Form("htemp41_%d",icent));
    for ( int i=0 ; i<=20 ; i++) {
      htemp41->SetBinContent(i,0);
      htemp41->SetBinError(i,0);
    }
    if (scaleByR) {
      htemp41->SetAxisRange(0,2,"Y");
      htemp41->SetYTitle("#frac{1}{N_{#gamma}} #frac{dN_{J#gamma}}{dx_{J#gamma}}");
    }
    htemp41->DrawCopy("hist");

    if (scaleByR)    {
      hxgj[kppdata13][icent+1]->Scale(hRpp2013[icent+1]->GetBinContent(1));
      std::cout <<"  Scaled by pp R :" << hRpp2013[icent+1]->GetBinContent(1) << std::endl;
      
      double x,y;
      rxhidata->GetPoint(icent, x,y);
      hxgj[khidata][icent]->Scale(y);
      std::cout << " scaled by PbPb R: " << y << std::endl;

      rxhimc->GetPoint(icent, x,y);
      hxgj[khimc][icent]->Scale(y);
      std::cout << " scaled by PbPb R: " << y << std::endl;
    }
    if(drawMC) hxgj[khimc][icent]->DrawCopy("hist same");

    if ( !mcOnly )    drawXSys(icent,hxgj[khidata][icent]);
    //   if ( icent==3){ }
    
    handsomeTH1(hxgj[kppdata13][icent+1]);
    hxgj[kppdata13][icent+1]->SetMarkerStyle(21);

    drawSys(hxgj[kppdata13][icent+1],ppSysX[icent],kGreen,3001);

    hxgj[kppdata13][icent+1]->Draw("same ");
    
    if ( !mcOnly )   hxgj[khidata][icent]->Draw("same");
    //onSun(0,0,2,0);
 
    // if ( icent == 2) {
    //   TLegend *leg0  = new TLegend(0.2796373,0.7545885,0.9742202,0.9937661,NULL,"brNDC");
    //     easyLeg(leg0,"");

CompareSpectra(Int_t part, Int_t charge, Int_t cent = -1, Int_t ratio = kFALSE, Int_t fitfunc = -1, Bool_t cutSpectra = kTRUE) 
{

  gROOT->LoadMacro("HistoUtils.C");

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

  LoadLibraries();
  AliBWFunc bwf;
  bwf.SetVarType(AliBWFunc::kdNdpt);
  TF1 *fFitFunc = NULL;

  switch (fitfunc) {
  case 0:
    gROOT->LoadMacro("SpectraAnalysis.C");
    fFitFunc = STAR_BlastWave("fBW", AliPID::ParticleMass(part), 0.9, 0.1, 1.);
    fBW->SetParLimits(3, 0.5, 1.5);
    //    fBW->FixParameter(3, 1.);
    break;
  case 1:
    fFitFunc = bwf.GetLevi(AliPID::ParticleMass(part), AliPID::ParticleMass(part), 5., 1000.);
    break;
  case 2:
    fFitFunc = bwf.GetBoltzmann(AliPID::ParticleMass(part), AliPID::ParticleMass(part), 100.);
    break;
  case 3:
    fFitFunc = bwf.GetMTExp(AliPID::ParticleMass(part),AliPID::ParticleMass(part) , 100.);
    break;
  case 4:
    fFitFunc = bwf.GetPTExp(AliPID::ParticleMass(part), 100.);
    break;
  case 5:
    fFitFunc = bwf.GetBGBW(AliPID::ParticleMass(part), 0.5, 0.1, 1.e6);
    break;
  case 6:
    fFitFunc = new TF1("fpol9", "pol9", 0., 5.0);
    break;
  }
  if (fFitFunc) fFitFunc->SetLineWidth(2);

  TFile *itssafile = TFile::Open(ratio ? itssaratiofilename : itssafilename);
  //  TFile *itstpcfile = TFile::Open(itstpcfilename);
  if (part / 3 == charge / 3)
    Char_t *tpctofratiofilename = tpctofratiofilenameA;
  else
    Char_t *tpctofratiofilename = tpctofratiofilenameB;
  TFile *tpctoffile = TFile::Open(ratio ? tpctofratiofilename : tpctoffilename);
  TFile *toffile = TFile::Open(ratio ? tofratiofilename : toffilename);
  //  TFile *hydrofile = TFile::Open(hydrofilename);

  TCanvas *cCanvas = new TCanvas("cCanvas");
  if (cent == -1) cCanvas->Divide(5, 2);
  TCanvas *cCanvasCombined = new TCanvas("cCanvasCombined");
  TCanvas *cCanvasRatio = new TCanvas("cCanvasRatio");
  if (cent == -1) cCanvasRatio->Divide(5, 2);
  TCanvas *cCanvasRatioComb = new TCanvas("cCanvasRatioComb");
  if (cent == -1) cCanvasRatioComb->Divide(5, 2);
  TCanvas *cCanvasRatioFit = new TCanvas("cCanvasRatioFit");
  if (cent == -1) cCanvasRatioFit->Divide(5, 2);
  TPad *curpad = NULL;
  TH1D *hITSsa, *hITSTPC, *hTPCTOF, *hTOF;
  TGraph *hHydro;
  TGraphErrors *gCombined[10];
  TProfile *pCombined[10];
  TH1D *hCombined[10];
  TH1D *hRatio_ITSsa_ITSTPC[10];
  TH1D *hRatio_ITSsa_TPCTOF[10];
  TH1D *hRatio_ITSTPC_TPCTOF[10];
  TH1D *hRatio_ITSTPC_TOF[10];
  TH1D *hRatio_TPCTOF_TOF[10];
  TH1D *hRatio_ITSsa_TOF[10];
  for (Int_t icent = 0; icent < 10; icent++) {
    if (cent != -1 && icent != cent) continue;
    gCombined[icent] = new TGraphErrors();
    pCombined[icent] = new TProfile(Form("pCombined_cent%d", icent), "", NptBins, ptBin);
    hCombined[icent] = new TH1D(Form("hCombined_cent%d", icent), "", NptBins, ptBin);
    TObjArray spectraArray;
    hITSsa = ratio ? GetITSsaRatio(itssafile, part, charge, icent, cutSpectra): GetITSsaSpectrum(itssafile, part, charge, icent, cutSpectra);
    //    hITSTPC = GetITSTPCSpectrum(itstpcfile, part, charge, icent, cutSpectra);
    hTPCTOF = ratio ? GetTPCTOFRatio(tpctoffile, part, charge, icent, cutSpectra) : GetTPCTOFSpectrum(tpctoffile, part, charge, icent, cutSpectra);
    hTOF = ratio ? GetTOFRatio(toffile, part, charge, icent, cutSpectra) : GetTOFSpectrum(toffile, part, charge, icent, cutSpectra);
    //    hHydro = GetHydroSpectrum(hydrofile, part, charge, icent);
    if (cent == -1)
      curpad = (TPad *)cCanvas->cd(icent + 1);
    else
      curpad = (TPad *)cCanvas->cd();
    if (!ratio)
      TH1D *hArea = new TH1D(Form("hArea_%d", icent), Form("%s (%s);p_{T} (GeV/c);#frac{d^{2}N}{dy dp_{t}};", partChargeName[part][charge], centName[icent]), 100, 0., 5.);
    else
      TH1D *hArea = new TH1D(Form("hArea_%d", icent), Form("%s (%s);p_{T} (GeV/c);#frac{d^{2}N}{dy dp_{t}};", "generic ratio", centName[icent]), 100, 0., 5.);
    hArea->Draw();
    Double_t minimum = 0.001;
    Double_t maximum = 1000.;
    if (hITSsa) {
      AddPointsToGraph(hITSsa, gCombined[icent]);
      AddPointsToProfile(hITSsa, pCombined[icent]);
      spectraArray.Add(hITSsa);
      hITSsa->DrawCopy("same");
      if (hITSsa->GetMaximum() > maximum) maximum = hITSsa->GetMaximum();
//.........这里部分代码省略.........

void drawDphi() { 
  
  TH1D* hdphi[5][5]; // [species][centrality] 
  TH1D* hxjg[5][5]; // [species][centrality] 
  TH1D* hxjgNorm[5][5];
  TFile* fff[5];
  
  fff[kHIDATA] = new TFile("ffFiles/photonTrackCorr_pbpbDATA_output_photonPtThr60_jetPtThr30_20130210_genLevel0.root");
  for ( int iSpecies =0 ; iSpecies<=5 ; iSpecies++) {
    if (iSpecies != kHIDATA)       continue;
    for ( int icent = 1 ; icent<=4 ; icent++) {
      hdphi[iSpecies][icent] = (TH1D*)fff[iSpecies]->Get(Form("jetDphi_icent%d_final",icent));
      hxjg[iSpecies][icent] = (TH1D*)fff[iSpecies]->Get(Form("jetXjg_icent%d_final",icent));
      hxjgNorm[iSpecies][icent] = (TH1D*)hxjg[iSpecies][icent]->Clone(Form("norm_%s",hxjg[iSpecies][icent]->GetName()));
      scaleInt(hxjgNorm[iSpecies][icent]);    
    }}
  
  fff[kHIMC]   = new TFile("ffFiles/photonTrackCorr_pbpbMC_output_photonPtThr60_jetPtThr30_20130210_genLevel0.root"); 
  for ( int iSpecies =0 ; iSpecies<=5 ; iSpecies++) {
    if (iSpecies != kHIMC)         continue;
    for ( int icent = 1 ; icent<=4 ; icent++) {
      hdphi[iSpecies][icent] = (TH1D*)fff[iSpecies]->Get(Form("jetDphi_icent%d_final",icent));
      hxjg[iSpecies][icent] = (TH1D*)fff[iSpecies]->Get(Form("jetXjg_icent%d_final",icent));
      hxjgNorm[iSpecies][icent] = (TH1D*)hxjg[iSpecies][icent]->Clone(Form("norm_%s",hxjg[iSpecies][icent]->GetName()));
      scaleInt(hxjgNorm[iSpecies][icent]);
    }}
  
  
  
  TCanvas* c1 = new TCanvas("c1","",1400,450);
  makeMultiPanelCanvas(c1,4,1,0.0,0.0,0.2,0.15,0.02);
  for ( int icent = 1 ; icent<=4 ; icent++) {
    c1->cd(5-icent);
    handsomeTH1(hdphi[kHIMC][icent],1);
    handsomeTH1(hdphi[kHIDATA][icent],2);
    hdphi[kHIMC][icent]->Scale(1./hdphi[kHIMC][icent]->Integral("width"));
    hdphi[kHIDATA][icent]->Scale(1./hdphi[kHIDATA][icent]->Integral("width"));
    hdphi[kHIMC][icent]->SetAxisRange(3.141592/3.,3.141592,"X");
    hdphi[kHIMC][icent]->SetAxisRange(0.01,5,"Y");
    hdphi[kHIMC][icent]->DrawCopy("hist");
    //  hdphi[kHIDATA][icent]->Draw("same");
    //    gPad->SetLogy();
    
    TF1 *fdphi = new TF1("fdphi","exp(-(TMath::Pi()-x)/[0]) / ([0]*(1-exp(-TMath::Pi()/[0]))) ",2.0*TMath::Pi()/3.0,TMath::Pi());
    fdphi->SetParName(0,"width");
    fdphi->SetParameter("width",0.22);
        
    float fitxmin=3.1415926*2./3;
    //    float fitxmin=2.5;
    hdphi[kHIMC][icent]->Fit("fdphi","0llm","",fitxmin,3.1415926);
    fdphi->SetLineWidth(1);
    fdphi->SetLineStyle(2);
    fdphi->SetLineColor(kBlue);
    float dphiWidth = fdphi->GetParameter("width");
    float dphiWidthErr = fdphi->GetParError(0);
    cout << " dphiWidth,dphiWidthErr = " << dphiWidth <<"   "<< dphiWidthErr << endl;
    fdphi->DrawCopy("same");
  }
  c1->SaveAs("pdfFiles/dPhi_MC_only.pdf");


  TCanvas* c2 = new TCanvas("c2","",1400,450);
  makeMultiPanelCanvas(c2,4,1,0.0,0.0,0.2,0.15,0.02);
  for ( int icent = 1 ; icent<=4 ; icent++) {
    c2->cd(5-icent);
    handsomeTH1(hxjgNorm[kHIMC][icent],1);
    handsomeTH1(hxjgNorm[kHIDATA][icent],2);
    hxjgNorm[kHIMC][icent]->SetAxisRange(-0.01,0.35,"Y");
    hxjgNorm[kHIMC][icent]->DrawCopy("");
    //    hxjgNorm[kHIDATA][5-icent]->Draw("same");
    int lowCent = centBin1[icent-1];    int highCent = centBin1[icent]-1;
    drawText(Form("%.0f%% - %.0f%%", float((float)lowCent*2.5), float((float)(highCent+1)*2.5)), 0.67,0.7,1);
  }
  TLegend *l12 = new TLegend(0.304355,0.755085,.8,0.9552542,NULL,"brNDC");
  easyLeg(l12,"p_{T}^{Jet}/p_{T}^{#gamma}");
  l12->AddEntry(hxjgNorm[kHIMC][1],"PYTHIA+HYDJET","pl");
  c2->cd(1);   l12->Draw();

  c2->SaveAs("pdfFiles/xgj_MC_only.pdf");

  TCanvas* c2MCover = new TCanvas("c2MCover","",500,500);
  bool drawn=false;
  for ( int icent = 1 ; icent<=4 ; icent++) {
    c2MCover->cd(icent);
    handsomeTH1(hxjgNorm[kHIMC][icent],ycolor[icent]);
    hxjgNorm[kHIMC][5-icent]->SetAxisRange(-0.01,0.38,"Y");
    if (drawn)  
      hxjgNorm[kHIMC][5-icent]->Draw("same");
    else       hxjgNorm[kHIMC][5-icent]->Draw("");
    
    
    drawn = true;
  }
  TLegend *l11 = new TLegend(0.1794355,0.7055085,1,0.9152542,NULL,"brNDC");
  easyLeg(l11,"p_{T}^{Jet}/p_{T}^{#gamma} of PYTHIA embedded in"); 
  l11->AddEntry(hxjgNorm[kHIMC][4],"HYDJET 50-100%","p");
  l11->AddEntry(hxjgNorm[kHIMC][3],"HYDJET 30-50%","p");
  l11->AddEntry(hxjgNorm[kHIMC][2],"HYDJET 10-30%","p");
  l11->AddEntry(hxjgNorm[kHIMC][1],"HYDJET 0 -10%","p");
  l11->Draw();
//.........这里部分代码省略.........

int FitSPE(TH1* spe, int ntriggers) 
{
  double fitmin, fitmax;
  double params[NPAR];
  //find the likely place for the single p.e. peak
  params[MEAN] = spe->GetMean();
  params[SIGMA] = spe->GetRMS();
  int bin = spe->FindBin(params[MEAN]);
  if(spe->GetMaximumBin() == 1){
    spe->Fit("gaus","Q0","",spe->GetBinCenter(bin), 
	     spe->GetBinCenter(bin+20));
  }
  else{
    spe->Fit("gaus","Q0","",spe->GetBinCenter(bin-10), 
	     spe->GetBinCenter(bin+10));
  }
  TF1* gaus = spe->GetFunction("gaus");
  if(gaus->GetParameter(1) > 0)
    params[MEAN] = gaus->GetParameter(1);
  if(gaus->GetParameter(2) > 0 && gaus->GetParameter(2) < params[SIGMA])
    params[SIGMA] = gaus->GetParameter(2);
  
  params[LAMBDA] = 0.5;
  //find the maximum fit range
  bin = spe->GetNbinsX();
  while(spe->GetBinContent(--bin) < 2) {}
  fitmax = spe->GetBinCenter(bin);
  
  spe->Fit("expo","Q0","",spe->GetBinCenter(bin), spe->GetBinCenter(bin-10));
  
  //find the best estimate for the exponentials
  //the first is from the first few bins
  bin = spe->FindBin(params[MEAN]) - 5;
  double compval = spe->GetBinContent(bin--);
  while(spe->GetBinContent(bin) < compval && spe->GetBinContent(bin) != 0){
    compval = spe->GetBinContent(bin--);
  }
  if(spe->GetBinContent(bin) == 0){
    fitmin = spe->GetBinCenter(bin+1);
  }
  else{
    ++bin;
    //now find the max of the turnover
    int minbin = bin;
    compval = spe->GetBinContent(bin--);
    while(spe->GetBinContent(bin) > compval && bin > 0 ){
      compval = spe->GetBinContent(bin--);
    }
    ++bin;
    spe->Fit("expo","Q0","",spe->GetBinCenter(bin), spe->GetBinCenter(minbin));
    params[AMP1] = spe->GetFunction("expo")->GetParameter(0);
    params[L1] = spe->GetFunction("expo")->GetParameter(1);
    fitmin = spe->GetBinCenter(bin+1);
  }
  
  //estimate the peak heights
  
  params[CONSTANT] = ntriggers * spe->GetBinWidth(1);
  fitmin = spe->GetBinCenter(4);
  
  TF1* spefunc = new TF1("spefunc", SPEFunc, fitmin, fitmax, NPAR);
  
  
  spefunc->SetParameters(params);
  for(int i=0; i<NPAR; i++)
    spefunc->SetParName(i, names[i]);
  spefunc->SetParLimits(CONSTANT, 0 , params[CONSTANT]*100.);
  spefunc->SetParLimits(LAMBDA, 0, 2);
  spefunc->SetParLimits(MEAN, std::max(params[MEAN]-params[SIGMA],0.), 
			params[MEAN]+params[SIGMA]);
  spefunc->SetParLimits(SIGMA, 0, spe->GetRMS());
    
  spefunc->SetParLimits(L1, -30, 0);
  spefunc->SetParLimits(AMP1, -30, 30);
  spefunc->SetParLimits(L2, -30, 0);
  spefunc->SetParLimits(AMP2, -30, 30);
    
  //spefunc->FixParameter(CONSTANT, ntriggers * spe->GetBinWidth(1));
  
  //return 0;
  spefunc->SetLineStyle(1);
  spefunc->SetLineColor(kBlue);
  int result = spe->Fit(spefunc,"MRNI");
  spefunc->DrawCopy("same");
  std::cout<<"Fit results: \n"<<
    "\t Chi2/NDF = "<<spefunc->GetChisquare()<<"/"<<spefunc->GetNDF()
	   <<"\n\t Prob = "<<spefunc->GetProb()<<std::endl;
  for(int i=0; i<NPAR; i++)
    params[i] = spefunc->GetParameter(i);
  TF1* background = new TF1("background",background_func,fitmin,fitmax,NPAR);
  background->SetLineColor(kRed);
  background->SetParameters(spefunc->GetParameters());
  background->DrawCopy("same");
  
  TF1* signal = new TF1("signal",signal_func,fitmin,fitmax,NPAR);
  signal->SetLineColor(kGreen); 
  signal->SetParameters(spefunc->GetParameters());
  signal->DrawCopy("same");
  
  TF1* apeak = new TF1("apeak","[0]*TMath::Gaus(x,[1],[2],1)",fitmin, fitmax);
//.........这里部分代码省略.........

//-------------------------------------------------------------------
void Alignment_1(int cluster_length_leftCut = 2,int cluster_length_rightCut = 7,float TotSignal_Cut = 5){
  gStyle->SetOptStat(11111111);
  gStyle->SetOptFit(111111111);
  //gAnaCombine->LinkBranches();
  ofstream outputlog("./AMS/AMSAlignment_1.txt",ios::out | ios::app);
  outputlog<<"\nfile: "<<gRootIOSvc->GetInputFileName()<<".\tCluster length: "<<cluster_length_leftCut<<"~"<<cluster_length_rightCut<<".\tTotal signale cut "<<TotSignal_Cut<<endl;

  TH1D *h_ladder_stripID[5][2];
  TH1D *h_Refer[5][2];
  char name[80];
  for(short lid =0;lid<5;++lid){
    for(short sid=0;sid<2;++sid){
      snprintf(name,80,"ladder_%d-side_%d",lid,sid);
      h_ladder_stripID[lid][sid] = new TH1D(name,name,500,0,1100);
      snprintf(name,80,"ladder_%d-side_%d VS ladder 0",lid,sid);
      h_Refer[lid][sid] = new TH1D(name,name,200,-200,200);
    }
    snprintf(name,80,"ladder_%d",lid);
    h_ladder_stripID[lid][0]->SetTitle(name);
    snprintf(name,80,"Seed Address ladder_%d",lid);
    h_Refer[lid][0]->SetTitle(name);
  }

  TChain *tree = gAnaCombine->GetTree();
  long nEvt = tree->GetEntries();
  std::vector<Cluster*> clusterIndex[5][2];
  for(long ievt=0;ievt<nEvt;++ievt){
    tree->GetEntry(ievt);
    short nCluster = gAnaCombine->fEvtAMS->GetEntriesFast();
    for(int iclu=0;iclu<nCluster;++iclu){
      Cluster *acluster = (Cluster*)gAnaCombine->fEvtAMS->At(iclu);
      if(acluster->GetTotSig()>TotSignal_Cut && cluster_length_leftCut<=acluster->length && acluster->length <= cluster_length_rightCut){
        int ladder = acluster->ladder;
        if(ladder==5){
          ladder = 2;
        }else if(ladder == 4){
          ladder = 3;
        }else if(ladder == 3){
          ladder = 4;
        }
        clusterIndex[ladder][acluster->side].push_back(acluster);
      }
    }
    if(clusterIndex[0][0].size() == 1 && clusterIndex[0][1].size() == 1){
      for(short il=0;il<5;++il){
        for(short is=0;is<2;++is){
          if(clusterIndex[il][is].size() == 1){
            h_ladder_stripID[il][is]->Fill(clusterIndex[il][is][0]->GetSeedAdd());
            h_Refer[il][is]->Fill(clusterIndex[il][is][0]->GetSeedAdd() - clusterIndex[0][is][0]->GetSeedAdd());
          }
          clusterIndex[il][is].clear();
        }
      }
    }
  }
  snprintf(name,80,"AMS Alignment_1_a: cluster length %d~%d, total signal cut %f",cluster_length_leftCut,cluster_length_rightCut,TotSignal_Cut);
  TCanvas *c1 = new TCanvas(name,name);
  c1->Divide(2,3);
  snprintf(name,80,"AMS Alignment_1_b: cluster length %d~%d, totalt signal cut %f",cluster_length_leftCut,cluster_length_rightCut,TotSignal_Cut);

  TF1 *gausFit = new TF1("GausFit","gaus",0,150);
  TCanvas *c2 = new TCanvas(name,name);
  c2->Divide(2,3);
  double mean=0.,sigma =0.;
  for(short lid=0;lid<5;++lid){
    c1->cd(lid+1);
    h_ladder_stripID[lid][0]->Draw();
    h_ladder_stripID[lid][1]->SetLineColor(6);
    h_ladder_stripID[lid][1]->Draw("same");

    c2->cd(lid+1);
    h_Refer[lid][0]->Draw();
    mean = h_Refer[lid][0]->GetMean(); sigma = h_Refer[lid][0]->GetRMS();
    gausFit->SetRange(mean-2.5*sigma,mean+2.5*sigma);
    h_Refer[lid][0]->Fit(gausFit,"R0Q");
    gausFit->DrawCopy("lsame");
    outputlog<<lid<<"  "<<gausFit->GetParameter(1)<<"  "<<gausFit->GetParameter(2)<<"    ";

    h_Refer[lid][1]->SetLineColor(6);
    h_Refer[lid][1]->Draw("same");
    mean = h_Refer[lid][1]->GetMean(); sigma = h_Refer[lid][1]->GetRMS();
    gausFit->SetRange(mean-2.5*sigma,mean+2.5*sigma);
    h_Refer[lid][1]->Fit(gausFit,"R0Q");
    gausFit->DrawCopy("lsame");
    outputlog<<gausFit->GetParameter(1)<<"  "<<gausFit->GetParameter(2)<<endl;
  }
}

//.........这里部分代码省略.........
				FitFuncGausSlices->SetParameters(hProfileY->GetBinContent(hProfileY->GetMaximumBin()),MaxValue,4);
				
				hProfileY->Fit(FitFuncGausSlices,"RNQ");
				
				FitFuncSlices->SetParameters(FitFuncGausSlices->GetParameter(0),FitFuncGausSlices->GetParameter(1),FitFuncGausSlices->GetParameter(2),10);
				
				FitFuncSlices->SetParLimits(0,0,10000);
				FitFuncSlices->SetParLimits(1,MaxValue-10,MaxValue+10);
				FitFuncSlices->SetParLimits(2,0,10);
				FitFuncSlices->SetParLimits(3,0,100);
				hProfileY->Fit(FitFuncSlices,"RNQ");
				
				cout <<MaxValue<<"  " << FitFuncSlices->GetParameter(1) << "   " << FitFuncSlices->GetParError(1) <<endl;
				hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually->SetBinContent(binX, FitFuncSlices->GetParameter(1));
				hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually->SetBinError(binX, FitFuncSlices->GetParError(1));
				
				
			}
			hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually->GetYaxis()->SetRangeUser(ChannelPeakPos-100,ChannelPeakPos+50);
			hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually->GetXaxis()->SetRangeUser(-0.05,0.9);
			TF1 * funcCorrConst = new TF1("funcCorrConst","pol0",-0.03,0.03);
			funcCorrConst->SetLineColor(kRed);
			TF1 * funcCorrPol = new TF1("funcCorrPol",PolPiecewise,-0.05,0.9,6);
			
			funcCorrPol->SetLineColor(kBlue);
			funcCorrPol->SetParameter(0,0.04);
			funcCorrPol->SetParameter(1,ChannelPeakPos);
			hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually->Fit(funcCorrConst,"R");
			hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually->Fit(funcCorrPol,"R+");
			TCanvas* can2=new TCanvas();
			gPad->SetLogz();
			
			hSpectrumTDeriMax1090Rel_EnergyChannel->Draw("colz");
			funcCorrPol->DrawCopy("same");
			for(int i = 7;i>0;i--)
				funcCorrPol->SetParameter(i,funcCorrPol->GetParameter(i)/funcCorrPol->GetParameter(1));
			funcCorrPol->Print();
			char buf[20];
			sprintf(buf, "funcCorrNorm_%d",ki);
			funcCorrPol->Write(buf,TObject::kOverwrite);
			
			funcCorrConst->Delete();
			funcCorrPol->Delete();
			
			
			//second correction via T1090 (after first)
			
			hSpectrumT1090_EnergyChannelCorr1->GetYaxis()->SetRangeUser(ChannelRangeMin,ChannelRangeMax);
			
			//TF1* FitFuncSlices = new TF1("FitFuncSlices","gaus(0)+pol0(3)",ChannelRangeMin,ChannelRangeMax);
			//FitFuncSlices->SetParameters(1000,ChannelPeakPos-10,4,10);
			//FitFuncSlices->SetParLimits(1,ChannelRangeMin,ChannelRangeMax);
			//FitFuncSlices->SetParLimits(2,0,10);
			//FitFuncSlices->SetParLimits(3,0,100);
	
			//gDirectory->ls();
			
			TH1D *hSpectrumT1090_EnergyChannelCorr1_MaxPosManually=new TH1D("hSpectrumT1090_EnergyChannelCorr1_MaxPosManually","",hSpectrumT1090_EnergyChannelCorr1->GetNbinsX(),0,1000);
			//cout <<hSpectrumTDeriMax1090_EnergyChannel_MaxPos->GetEntries()<< endl;
			for(int binX = hSpectrumT1090_EnergyChannelCorr1->GetXaxis()->FindBin(30);binX <= hSpectrumT1090_EnergyChannelCorr1->GetXaxis()->FindBin(500);binX++)
			{
				cout << "binx " << binX << endl;
				TH1D *hProfileY =hSpectrumT1090_EnergyChannelCorr1->ProjectionY("_py",binX,binX);
				//hProfileY->Draw();
				//return 0;
				double MaxValue=hProfileY->GetBinCenter(hProfileY->GetMaximumBin());

//.........这里部分代码省略.........
				// Find next regular term only
				bool foundRegTerm = false;
				for (int jTerm = iTerm + 1; jTerm < nTerms && !foundRegTerm; ++jTerm) {
					const Student::Enrollment jEnrollment = student->Enrollments()[jTerm];
					if (!MyFunctions::regularSemester(jEnrollment.term)) continue;
					foundRegTerm = true;
					for (Student::Grade jGrade : jEnrollment.grades) {
						if (jGrade.course == iGrade.course) continue;
						if (!MyFunctions::ValidGrade(jGrade.grade)) continue;
						
						double prediction_j = student->CourseGradePrediction(jGrade, Student::DISTRIBUTION);
						double delta_j = jGrade.quality - prediction_j;
						
						corrMap[std::make_pair(iGrade.course, jGrade.course)].Add(delta_i, delta_j);
						++nPairAll;
						
					}
				}
			} 
		}
		
	}
	myBenchmark->Stop("Main Loop");
	
	std::cout << "nPairAll     = " << nPairAll << std::endl;
	std::cout << "Unique Pairs = " << corrMap.size() << std::endl;
	
	TH1D* rHist = new TH1D("rHist", "Correlation Coefficient, #rho", 120, -1.2, 1.2);
	TH1D* pHist = new TH1D("pHist", "Probablity Distribution", 100, 0., 1.);
	TH1D* nHist = new TH1D("nHist", "Number of entries", 100, 0., 2000.);
	TH2D* pVrHist = new TH2D("pVrHist", "Prob vs. #rho", 100, -1., 1., 100, 0., 1.);
	
	myBenchmark->Start("Prune");
	for (auto iter = corrMap.begin(); iter != corrMap.end();) {
		if (iter->second.n() < nCut) {
			corrMap.erase(iter++);
			continue;
		}
		double p = iter->second.p();
		double r = iter->second.r();

		// Test for nan?
		if (p != p) {
			std::cout << "Found p = nan: n = " << iter->second.n() << std::endl;
			corrMap.erase(iter++);
			continue;
		}
		if (p < 0.) {
			corrMap.erase(iter++);
			continue;
		}
		rHist->Fill(r);
		pHist->Fill(p);
		nHist->Fill(iter->second.n());
		pVrHist->Fill(r, p);
		if (p < 1. - prob && p > prob) {
			corrMap.erase(iter++);
		}
		else {
			std::cout << "r = " << iter->second.r() << ", p = " << p << std::endl;
			++iter;
		}
	}
	myBenchmark->Stop("Prune");
	std::cout << "Post Cut     = " << corrMap.size() << std::endl;
	
	myBenchmark->Start("Sort");
	std::vector<std::pair<std::pair<TString, TString>, CorrelationCalculator>> corrVec(corrMap.begin(), corrMap.end());
	std::sort(corrVec.begin(), corrVec.end(), &sortFunc);
	myBenchmark->Stop("Sort");
	
	int printTop = 50;
	int printed = 0;
	for (auto const& entry : corrVec) {
		std::cout << entry.first.first << " : " << entry.first.second << "\t, n = " << entry.second.n() << "\t, r = " << entry.second.r() 
			<< "\t, p = " << entry.second.p() << std::endl;
		++printed;
		if (printed >= printTop) break;
	}

	TCanvas* c1 = new TCanvas("c1", "Grade Correlation", 1600, 1200);
	c1->Divide(2,2);
	c1->cd(1);
	TF1* myGaus = new TF1("myGaus", "gaus", -1., 1.);
	myGaus->SetParameters(600., 0., 0.2);
	myGaus->FixParameter(1, 0.);
	rHist->Fit(myGaus, "0B", "", -1., 0.);
	rHist->DrawCopy();
	myGaus->DrawCopy("SAME");
	c1->cd(2);
	pHist->DrawCopy();
	c1->cd(3);
	nHist->DrawCopy();
	c1->cd(4);
	pVrHist->DrawCopy();
	
	float rt, cp;
	myBenchmark->Summary(rt, cp);
	
}