C++ TH2D::FindBin方法代码示例

void runcorr(int filenum)
{


  string buffer;
  vector<string> listoffiles;
  int nlines = 0;
  ifstream infile("/net/hisrv0001/home/dav2105/corrana/makecorrhists/oldstuff/sortedforests.txt");

  if (!infile.is_open()) {
    cout << "Error opening file. Exiting." << endl;
    return;
  } else {
    while (!infile.eof()) {
      infile >> buffer;
      listoffiles.push_back(buffer);
      nlines++;
    }
  }


  corrana(listoffiles[filenum].data());
  
  double pttriglow[] = {1,4};
  double pttrighigh[] = {2,8};
  double ptasslow[] = {1,2};
  double ptasshigh[] = {2,3};
  int centmin[] = {0,4,8,12,16,20,24,28,32};
  int centmax[] = {40,8,12,16,20,24,28,32,36};
  TFile * outf = new TFile(Form("corrhists_%d.root",filenum),"recreate");
  
  // for(int i = 0 ; i < 9 ; ++i)
  // for(int i = 0 ; i < 2 ; ++i)
  for(int i = 0 ; i < 1 ; ++i)
  {
    cout<<"pt iteration "<<i<<endl;
    // for(int cent = 0 ; cent < 3 ; ++cent)
    // for(int cent = 0 ; cent < 9 ; ++cent)
    for(int cent = 0 ; cent < 1 ; ++cent)
    {
      cout<<"cent iteration "<<cent<<endl;
      TH2D * ttsig = TrackTrackSignal(pttriglow[i],pttrighigh[i],ptasslow[i],ptasshigh[i],centmin[cent],centmax[cent]);
      TH2D * ttbak = TrackTrackBackground(pttriglow[i],pttrighigh[i],ptasslow[i],ptasshigh[i],centmin[cent],centmax[cent]);
      // TCanvas * c1 = new TCanvas();
      // ttsig->Draw("surf1");
      // TCanvas * c2 = new TCanvas();
      // ttbak->Draw("surf1");
      // TCanvas * c3 = new TCanvas();
      TH2D * ttcorr = (TH2D*)ttsig->Clone(Form("corr_trg%d_%d_ass%d_%d_cmin%d_cmax%d",(int)pttriglow[i],(int)pttrighigh[i],(int)ptasslow[i],(int)ptasshigh[i],centmin[cent],centmax[cent]));
      ttcorr->Divide(ttbak);
      ttcorr->Scale(ttbak->GetBinContent(ttbak->FindBin(0,0)));
      ttcorr->GetXaxis()->SetRange(ttcorr->GetXaxis()->FindBin(-4.0),ttcorr->GetXaxis()->FindBin(4.0));
      ttcorr->GetYaxis()->SetRange(ttcorr->GetYaxis()->FindBin(-3.1415926/2.0),ttcorr->GetYaxis()->FindBin(3*3.1415926/2.0));
      // ttcorr->Draw("surf1");
    }
  }
  
  outf->Write();
  outf->Close();
}

void mcruncorr(int filenum = 0)
{


  string buffer;
  vector<string> listoffiles;
  int nlines = 0;
  ifstream infile("/net/hisrv0001/home/dav2105/corrana/makecorrhists/franksorted.txt");

  if (!infile.is_open()) {
    cout << "Error opening file. Exiting." << endl;
    return;
  } else {
    while (!infile.eof()) {
      infile >> buffer;
      listoffiles.push_back(buffer);
      nlines++;
    }
  }

  // int filenum = 13;
  
  cout<<"opening: "<<listoffiles[filenum].data()<<endl;
  // alicecorrana("/mnt/hadoop/cms/store/user/velicanu/mergedv1_sortedforest/mergesortv1_2.root");
  mccorrana(listoffiles[filenum].data());
  
  // double leadingjetptlow[] =     {0  ,70 ,70 ,120,90 ,60 ,120,90 ,60 };
  // double leadingjetpthigh[] =    {300,90 ,300,300,300,300,300,300,300};
  // double subleadingjetptlow[] =  {70 ,0  ,50 ,50 ,50 ,50 ,50 ,50 ,50 };
  // double subleadingjetpthigh[] = {90 ,300,300,300,300,300,300,300,300};
  // double ptasslow[] = {1,1,1,2,2,2,3,3,3};
  // double ptasshigh[] = {2,2,2,3,3,3,4,4,4};
  // int centmin[] = {0,4,8,16,24};
  // int centmax[] = {4,8,16,24,28};
  
  
  
  double leadingjetptlow[] =     {0  ,70 ,0  ,70 ,0  ,70 ,90 ,90 ,90 };
  double leadingjetpthigh[] =    {300,90 ,300,90 ,300,90 ,300,300,300};
  double subleadingjetptlow[] =  {70 ,0  ,70 ,0  ,70 ,0  ,50 ,50 ,50 };
  double subleadingjetpthigh[] = {90 ,300,300,300,300,300,300,300,300};
  double ptasslow[] =            {1  ,1  ,2  ,2  ,3  ,3  ,1  ,2  ,3  };
  double ptasshigh[] =           {2  ,2  ,3  ,3  ,4  ,4  ,2  ,3  ,4  };
  int centmin[] = {0,4,8,12,16,20};
  int centmax[] = {4,8,12,16,20,24};
  TFile * outf = new TFile(Form("frank_pf3_%d.root",filenum),"recreate");
  
  for(int i = 6 ; i < 9 ; ++i)
  {
    cout<<"pt iteration "<<i<<endl;
    // for(int cent = 0 ; cent < 3 ; ++cent)
    // for(int cent = 0 ; cent < 9 ; ++cent)
    for(int cent = 0 ; cent < 6 ; ++cent)
    {
      cout<<"cent iteration "<<cent<<endl;
      TH2D * ljtsig = JetTrackSignal(0, leadingjetptlow[i], leadingjetpthigh[i] , subleadingjetptlow[i] , subleadingjetpthigh[i] , ptasslow[i] , ptasshigh[i], centmin[cent], centmax[cent]);
      TH2D * ljtbak = JetTrackBackground(0, leadingjetptlow[i], leadingjetpthigh[i] , subleadingjetptlow[i] , subleadingjetpthigh[i] , ptasslow[i] , ptasshigh[i], centmin[cent], centmax[cent]);
      TH2D * ljtcorr = (TH2D*)ljtsig->Clone(Form("corr_leadingjet%d_%d_ass%d_%d_cmin%d_cmax%d",(int)leadingjetptlow[i],(int)leadingjetpthigh[i],(int)ptasslow[i],(int)ptasshigh[i],centmin[cent],centmax[cent]));
      ljtcorr->Divide(ljtbak);
      ljtcorr->Scale(ljtbak->GetBinContent(ljtbak->FindBin(0,0)));
      ljtcorr->GetXaxis()->SetRange(ljtcorr->GetXaxis()->FindBin(-1.6),ljtcorr->GetXaxis()->FindBin(1.6));
      ljtcorr->GetYaxis()->SetRange(ljtcorr->GetYaxis()->FindBin(-3.1415926/2.0),ljtcorr->GetYaxis()->FindBin(3*3.1415926/2.0));
      
        TH2D * sljtsig = JetTrackSignal(1, leadingjetptlow[i], leadingjetpthigh[i] , subleadingjetptlow[i] , subleadingjetpthigh[i] , ptasslow[i] , ptasshigh[i], centmin[cent], centmax[cent]);
      TH2D * sljtbak = JetTrackBackground(1, leadingjetptlow[i], leadingjetpthigh[i] , subleadingjetptlow[i] , subleadingjetpthigh[i] , ptasslow[i] , ptasshigh[i], centmin[cent], centmax[cent]);
      TH2D * sljtcorr = (TH2D*)sljtsig->Clone(Form("corr_subleadingjet%d_%d_ass%d_%d_cmin%d_cmax%d",(int)leadingjetptlow[i],(int)leadingjetpthigh[i],(int)ptasslow[i],(int)ptasshigh[i],centmin[cent],centmax[cent]));
      sljtcorr->Divide(sljtbak);
      sljtcorr->Scale(sljtbak->GetBinContent(sljtbak->FindBin(0,0)));
      sljtcorr->GetXaxis()->SetRange(sljtcorr->GetXaxis()->FindBin(-1.6),sljtcorr->GetXaxis()->FindBin(1.6));
      sljtcorr->GetYaxis()->SetRange(sljtcorr->GetYaxis()->FindBin(-3.1415926/2.0),sljtcorr->GetYaxis()->FindBin(3*3.1415926/2.0));
      
      
    }
  }
  outf->Write();
  outf->Close();
}

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

            ring = (CbmRichRing*)fRichRings->At(iRing);
            if (!ring) continue;
            match = (CbmRichRingMatch*)fRichMatches->At(iRing);
            if (!match) continue;

            Int_t trackId = match->GetMCTrackID();
            if (trackId == -1) continue;
            if (trackId > fMCTracks->GetEntries()) continue;

            CbmMCTrack* mcTrack = (CbmMCTrack*)fMCTracks->At(trackId);
            if (!mcTrack) continue;
            Int_t pdg = TMath::Abs(mcTrack->GetPdgCode());
            Int_t motherId = mcTrack->GetMotherId();
            if (pdg != 11) continue;
            if (motherId != -1) continue;

            Double_t axisA = ring->GetAaxis();
            Double_t axisB = ring->GetBaxis();
            if (axisA > fMaxAaxis || axisB > fMaxAaxis) continue;
            if (axisA < fMinAaxis || axisB < fMinAaxis) continue;

            Double_t radius = ring->GetRadius();
            Double_t centerX = ring->GetCenterX();
            Double_t centerY = ring->GetCenterY();
            Double_t axisAbefore = ring->GetAaxis();
            Double_t axisBbefore = ring->GetBaxis();
            fh_Abefore->Fill(axisAbefore);
            fh_Bbefore->Fill(axisBbefore);

            Double_t axisA = ring->GetAaxis();
            Double_t axisB = ring->GetBaxis() ;

            axisA += mapaxisAXY->GetBinContent(mapaxisAXY->FindBin(centerX,centerY));
            axisB += mapaxisBXY->GetBinContent(mapaxisBXY->FindBin(centerX,centerY));

            fh_A->Fill(axisA);
            fh_B->Fill(axisB);
        } //iRing
    }//iEvent


  //  gStyle->SetOptStat(0);
    c3 = new TCanvas("c3","c3",10,10,600,600);
    c3->Divide(2,2);
    c3->cd(1);
    fh_Abefore->Scale(1./fh_Abefore->Integral());
    fh_Abefore->SetMaximum(fh_Abefore->GetMaximum()*1.3);
    fh_Abefore->Draw();
    fh_Abefore->SetAxisRange(fMinAaxis, fMaxAaxis);
    fh_Abefore->Fit("gaus");
    Double_t sigmaAb = fh_Abefore->GetFunction("gaus")->GetParameter("Sigma");
    char sigmaTxtAb[30];
    sprintf(sigmaTxtAb,"sigma = %.3f",sigmaAb);
    TText* txtAb = new TText(4.3, fh_Abefore->GetMaximum()*0.85, sigmaTxtAb);
    txtAb->SetTextSize(0.1);
    txtAb->Draw();
    gPad->SetGridx(true);
    gPad->SetGridy(true);

    c3->cd(2);
    fh_Bbefore->Scale(1./fh_Bbefore->Integral());
    fh_Bbefore->SetMaximum(fh_Bbefore->GetMaximum()*1.3);
    fh_Bbefore->Draw();
    fh_Bbefore->SetAxisRange(fMinAaxis, fMaxAaxis);
    fh_Bbefore->Fit("gaus");

void residualAlignment(TH2D* residualX, TH2D* residualY, double& offsetX,
                       double& offsetY, double& rotation,
                       double relaxation, bool display)
{
  assert(residualX && residualY && "Processors: can't perform residual alignment without histograms");

  rotation = 0;
  offsetX = 0;
  offsetY = 0;
  double angleWeights = 0;
  double fitChi2 = 0;

  for (int axis = 0; axis < 2; axis++)
  {
    TH2D* hist = 0;
    if (axis) hist = residualX;
    else      hist = residualY;

    // Project the histogram and fit with a gaussian to center the sensor
    TH1D* project = hist->ProjectionX("ResidualProjetion", 1, hist->GetNbinsY());
    project->SetDirectory(0);

    double sigma = project->GetBinWidth(1);
    double mean = 0;
    fitGaussian(project, mean, sigma, false);

    if (axis) offsetX = mean;
    else      offsetY = mean;

    delete project;

    std::vector<double> ptsX;
    std::vector<double> ptsY;
    std::vector<double> ptsErr;

    const unsigned int numSlices = hist->GetNbinsY();

    for (Int_t row = 1; row <= (int)numSlices; row++)
    {
      TH1D* slice = hist->ProjectionX("ResidualSlice", row, row);
      slice->SetDirectory(0);

      double mean = 0;
      double sigma = 0;
      double factor = 0;
      double background = 0;

      if (slice->Integral() < 1) { delete slice; continue; }
      fitGaussian(slice, mean, sigma, factor, background, false);

      const double sliceMin = slice->GetBinCenter(1);
      const double sliceMax = slice->GetBinCenter(slice->GetNbinsX());
      delete slice;

      // Quality assurance

      // Sigma is contained in the slice's range
      if (sigma > (sliceMax - sliceMin)) continue;
      // Mean is contained in the slice's range
      if (mean > sliceMax || mean < sliceMin) continue;
      // Peak is contains sufficient events
      if (factor < 100) continue;
      // Sufficient signal to noise ratio
      if (factor / background < 10) continue;

      // Get the total number of events in the gaussian 1 sigma
      Int_t sigRangeLow = hist->FindBin(mean - sigma);
      Int_t sigRangeHigh = hist->FindBin(mean + sigma);

      double sigRangeTotal = 0;
      for (Int_t bin = sigRangeLow; bin <= sigRangeHigh; bin++)
        sigRangeTotal += hist->GetBinContent(bin);

      // 2 * 1 sigma integral shoudl give ~ area under gaussian
      sigma /= sqrt(2 * sigRangeTotal);

      ptsX.push_back(hist->GetYaxis()->GetBinCenter(row));
      ptsY.push_back(mean);
      ptsErr.push_back(sigma);
    }

    if (ptsX.size() < 3) continue;

    std::vector<double> yvals = ptsY;
    std::sort(yvals.begin(), yvals.end());
    const double median = yvals[yvals.size()/2];
    double avgDeviation = 0;
    for (unsigned int i = 0; i < yvals.size(); i++)
      avgDeviation += fabs(yvals[i] - median);
    avgDeviation /= (double)yvals.size();

    std::vector<double> ptsXGood;
    std::vector<double> ptsYGood;
    std::vector<double> ptsErrGood;

    for (unsigned int i = 0; i < ptsX.size(); i++)
    {
      if (fabs(ptsY[i] - median) > 1.5*avgDeviation) continue;
      ptsXGood.push_back(ptsX[i]);
      ptsYGood.push_back(ptsY[i]);
//.........这里部分代码省略.........

void runcorr(int condor_iter, int trackqual)
{

  const int nptbins = 2;
  const int ncentbins = 1;
  const int najbins = 1;
  
  string buffer;
  vector<string> listoffiles;
  int nlines = 0;
  // ifstream infile("sortedforests.txt");
  // ifstream infile("doeproposalforests.txt");
  // ifstream infile("ppb_unmerged_minbias_forests.txt");
  ifstream infile("HIRun2013-PromptReco-v1-HLT_PAPixelTracks_Multiplicity190_v1-forest-v2.txt");

  if (!infile.is_open()) {
    cout << "Error opening file. Exiting." << endl;
    return;
  } else {
    while (!infile.eof()) {
      infile >> buffer;
      listoffiles.push_back(buffer);
      nlines++;
    }
  }
  bool dostdhists = (condor_iter%(nptbins * ncentbins * najbins) == 0);
  int ptbin = condor_iter % nptbins;
  int centbin = (condor_iter / nptbins) % ncentbins;
  int ajbin = (condor_iter / (nptbins * ncentbins)) % najbins;
  int filenum = (condor_iter / (nptbins * ncentbins * najbins));
  cout << "ipt: " << ptbin << " icent: " << centbin << " iaj: " << ajbin << " filenum: " << filenum << " dostdhists: " << dostdhists << " condor_iter "<<condor_iter<< endl;
 

  int nmin = 190 , nmax = 1000;
  //! for first iteration of forest production
  /*
  if(filenum==0) { nmin = 110 ; nmax = 1000; }
  if(filenum==1) { nmin = 90  ; nmax = 110 ; }
  if(filenum>1 ) { nmin = 0   ; nmax = 35  ; }
  if(filenum>9 ) { nmin = 35  ; nmax = 90  ; }
  */
  //! for second iteration of forest production
  // /*
  // if(filenum<24 ) { nmin = 90  ; nmax = 110 ; }
  // if(filenum<22 ) { nmin = 35  ; nmax = 90  ; }
  // if(filenum<12 ) { nmin = 110 ; nmax = 1000; }
  // if(filenum<10 ) { nmin = 0   ; nmax = 35  ; }
  // */
  //! for second iteration of forest production
  /*
  if(filenum<26 ) { nmin = 90  ; nmax = 110 ; }
  if(filenum<23 ) { nmin = 35  ; nmax = 90  ; }
  if(filenum<13)  { nmin = 110 ; nmax = 1000; }
  if(filenum<10)  { nmin = 0   ; nmax = 35  ; }
  */
  corrana(listoffiles[filenum].data(),trackqual);
  
  
  double pttriglow[] =  {1    ,1    ,1 ,3, 1};
  double pttrighigh[] = {3    ,3    ,3 ,4, 3};
  double ptasslow[] =   {0.25 ,2 ,0.25 ,3, 1};
  double ptasshigh[] =  {0.5  ,3    ,3 ,4, 3};
  
  
  int centmin[] = {1,0,4,8,12,16,20,24,28,32};
  int centmax[] = {1,41,8,12,16,20,24,28,32,36};
  TFile * outf = new TFile(Form("corrhists_trkhfminus_trkqaul%d_nmin%d_nmax%d_ptmin%d_ptmax%d_hfmin%d_hfmax%d_%d.root",trackqual,nmin,nmax,(int)pttriglow[ptbin],(int)pttrighigh[ptbin],(int)ptasslow[ptbin],(int)ptasshigh[ptbin],filenum),"recreate");
  
  int i = 0;
  int cent = 0;

  cout<<"cent iteration "<<cent<<endl;
  TH2D * ttsig = HFTrackSignal(pttriglow[ptbin],pttrighigh[ptbin],ptasslow[ptbin],ptasshigh[ptbin],centmin[cent],centmax[cent],nmin,nmax);
  TH2D * ttbak = HFTrackBackground(pttriglow[ptbin],pttrighigh[ptbin],ptasslow[ptbin],ptasshigh[ptbin],centmin[cent],centmax[cent],nmin,nmax);

  // TH1I * hntottrig = new TH1I(Form("nttottrig_trg%d_%d_ass%d_%d_cmin%d_cmax%d",(int)pttriglow[ptbin],(int)pttrighigh[ptbin],(int)ptasslow[ptbin],(int)ptasshigh[ptbin],centmin[cent],centmax[cent]),"",1,0.5,1.5);
  // int myntottrig = GetNTotTrig();
  // hntottrig->Fill(1,myntottrig);
  // cout<<"ntottrig: "<<myntottrig<<endl;

  TH2D * ttcorr = (TH2D*)ttsig->Clone(Form("corr_trg%d_%d_ass%d_%d_cmin%d_cmax%d",(int)pttriglow[ptbin],(int)pttrighigh[ptbin],(int)ptasslow[ptbin],(int)ptasshigh[ptbin],centmin[cent],centmax[cent]));
  ttcorr->Divide(ttbak);
  ttcorr->Scale(ttbak->GetBinContent(ttbak->FindBin(4,0)));
  ttcorr->Scale(1/0.0594998609); //! bin width
  ttcorr->GetXaxis()->SetRange(ttcorr->GetXaxis()->FindBin(1),ttcorr->GetXaxis()->FindBin(7));
  ttcorr->GetYaxis()->SetRange(ttcorr->GetYaxis()->FindBin(-3.1415926/2.0),ttcorr->GetYaxis()->FindBin(3*3.1415926/2.0));
  
  outf->Write();
  outf->Close();
}

void makePFJetResolutionMCtoDataCorrLUT()
{
  Double_t absEtaBins[] = { -0.1, 1.1, 1.7, 2.3, 9.9 };
  const Int_t numAbsEtaBins = 4;
  
  // numbers from JME-10-014 PAS
  Double_t mcToDataCorr_dijet2010[]              = { 1.04,  1.02,  0.91,  1.04  };
  Double_t mcToDataCorrStatErr_dijet2010[]       = { 0.012, 0.029, 0.044, 0.060 };
  Double_t mcToDataCorrSysErrUp_dijet2010[]      = { 0.052, 0.070, 0.085, 0.161 };
  Double_t mcToDataCorrSysErrDown_dijet2010[]    = { 0.049, 0.068, 0.081, 0.153 };

  Double_t mcToDataCorr_gammajet2010[]           = { 1.07,  1.10,  1.07,  1.18  };
  Double_t mcToDataCorrStatErr_gammajet2010[]    = { 0.020, 0.031, 0.048, 0.062 };
  Double_t mcToDataCorrSysErrUp_gammajet2010[]   = { 0.014, 0.009, 0.035, 0.017 };
  Double_t mcToDataCorrSysErrDown_gammajet2010[] = { 0.026, 0.023, 0.029, 0.057 };

  // numbers from AN-11/330
  Double_t mcToDataCorr_dijet2011[]              = { 1.055, 1.096, 1.134, 1.288 };
  Double_t mcToDataCorrStatErr_dijet2011[]       = { 0.012, 0.017, 0.035, 0.127 };
  Double_t mcToDataCorrSysErrUp_dijet2011[]      = { 0.061, 0.068, 0.097, 0.183 };
  Double_t mcToDataCorrSysErrDown_dijet2011[]    = { 0.060, 0.067, 0.095, 0.181 };

  TH2D* pfJetResolutionMCtoDataCorrLUT = new TH2D("pfJetResolutionMCtoDataCorrLUT", 
						  "pfJetResolutionMCtoDataCorrLUT", numAbsEtaBins, absEtaBins, 1, -0.01, 1.e+6);

  for ( Int_t iAbsEtaBin = 0; iAbsEtaBin < numAbsEtaBins; ++iAbsEtaBin ) {
    Double_t mcToDataCorrSysErr_dijet2010 = ( mcToDataCorr_dijet2010[numAbsEtaBins] < mcToDataCorr_gammajet2010[numAbsEtaBins] ) ?
      mcToDataCorrSysErrUp_dijet2010[iAbsEtaBin] : mcToDataCorrSysErrDown_dijet2010[iAbsEtaBin];
    Double_t mcToDataCorrTotErr2_dijet2010 = 
      square(mcToDataCorrStatErr_dijet2010[iAbsEtaBin]) + square(mcToDataCorrSysErr_dijet2010);
    Double_t weight_dijet2010 = 1./mcToDataCorrTotErr2_dijet2010;

    Double_t mcToDataCorrSysErr_gammajet2010 = ( mcToDataCorr_gammajet2010[numAbsEtaBins] < mcToDataCorr_dijet2010[numAbsEtaBins] ) ?
      mcToDataCorrSysErrUp_gammajet2010[iAbsEtaBin] : mcToDataCorrSysErrDown_gammajet2010[iAbsEtaBin];
    Double_t mcToDataCorrTotErr2_gammajet2010 = 
      square(mcToDataCorrStatErr_gammajet2010[iAbsEtaBin]) + square(mcToDataCorrSysErr_gammajet2010);
    Double_t weight_gammajet2010 = 1./mcToDataCorrTotErr2_gammajet2010;

    Double_t mcToDataCorr_2010 = 
      (weight_dijet2010*mcToDataCorr_dijet2010[iAbsEtaBin] + weight_gammajet2010*mcToDataCorr_gammajet2010[iAbsEtaBin])/
      (weight_dijet2010 + weight_gammajet2010);
    Double_t mcToDataCorrTotErr_2010 = TMath::Sqrt(1./(weight_dijet2010 + weight_gammajet2010));
    std::cout << "eta = " << absEtaBins[iAbsEtaBin] << "-" << absEtaBins[iAbsEtaBin + 1] << ":"
	      << " MC/Data correction = " << mcToDataCorr_2010 << " +/- " << mcToDataCorrTotErr_2010 << std::endl;

    Int_t binIndex = pfJetResolutionMCtoDataCorrLUT->FindBin(0.5*(absEtaBins[iAbsEtaBin] + absEtaBins[iAbsEtaBin + 1]), 1.);
    pfJetResolutionMCtoDataCorrLUT->SetBinContent(binIndex, mcToDataCorr_2010);
    pfJetResolutionMCtoDataCorrLUT->SetBinError(binIndex, mcToDataCorrTotErr_2010);    
  }
/*
  for ( Int_t iAbsEtaBin = 0; iAbsEtaBin < numAbsEtaBins; ++iAbsEtaBin ) {
    Double_t mcToDataCorrSysErr_dijet2010 = ( mcToDataCorr_dijet2010[numAbsEtaBins] < mcToDataCorr_dijet2011[numAbsEtaBins] ) ?
      mcToDataCorrSysErrUp_dijet[iAbsEtaBin] : mcToDataCorrSysErrDown_dijet2010[iAbsEtaBin];
    Double_t mcToDataCorrTotErr2_dijet2010 = 
      square(mcToDataCorrStatErr_dijet2010[iAbsEtaBin]) + square(mcToDataCorrSysErr_dijet2010);

    Double_t mcToDataCorrSysErr_dijet2011 = ( mcToDataCorr_dijet2011[numAbsEtaBins] < mcToDataCorr_dijet2010[numAbsEtaBins] ) ?
      mcToDataCorrSysErrUp_dijet2011[iAbsEtaBin] : mcToDataCorrSysErrDown_dijet2011[iAbsEtaBin];
    Double_t mcToDataCorrTotErr2_dijet2011 = 
      square(mcToDataCorrStatErr_dijet2011[iAbsEtaBin]) + square(mcToDataCorrSysErr_dijet2011);

    Double_t ratio_dijet2011to2010 = (mcToDataCorr_dijet2011[iAbsEtaBin]/mcToDataCorr_dijet2010[iAbsEtaBin]);
    Double_t ratioErr_dijet2011to2010 = 
      ratio_dijet2011to2010*TMath::Sqrt(mcToDataCorrTotErr2_dijet2010/square(mcToDataCorr_dijet2011[iAbsEtaBin])
				      + mcToDataCorrTotErr2_dijet2011/square(mcToDataCorr_dijet2010[iAbsEtaBin]));
    std::cout << "eta = " << absEtaBins[iAbsEtaBin] << "-" << absEtaBins[iAbsEtaBin + 1] << ":"
	      << " dijet2011/dijet2010 correction = " << ratio_dijet2011to2010
	      << " +/- " << ratioErr_dijet2011to2010 << std::endl;
  }
 */
/*
  TFile* outputFile = new TFile("../data/pfJetResolutionMCtoDataCorrLUT.root", "RECREATE");
  pfJetResolutionMCtoDataCorrLUT->Write();
  delete outputFile;
 */
}

void stdruncorr(int filenum = 0)
{


  string buffer;
  vector<string> listoffiles;
  int nlines = 0;
  ifstream infile("/net/hisrv0001/home/dav2105/corrana/makecorrhists/mcsorted.txt");

  if (!infile.is_open()) {
    cout << "Error opening file. Exiting." << endl;
    return;
  } else {
    while (!infile.eof()) {
      infile >> buffer;
      listoffiles.push_back(buffer);
      nlines++;
    }
  }

  
  cout<<"opening: "<<listoffiles[filenum].data()<<endl;
  stdcorrana(listoffiles[filenum].data());
  double leadingjetptlow[] =     {100,100,100,100,100,100,100,100,100};
  double leadingjetpthigh[] =    {120,120,120,300,300,300,300,300,300};
  double subleadingjetptlow[] =  {50 ,50 ,50 ,50 ,50 ,50 ,50 ,50 ,50 };
  double subleadingjetpthigh[] = {120,120,120,300,300,300,300,300,300};
  double ptasslow[] =            {2  ,3  ,5  ,2  ,3  ,3  ,1  ,2  ,3  };
  double ptasshigh[] =           {3  ,5  ,8  ,3  ,4  ,4  ,2  ,3  ,4  };


  int centmin[] = {0,4,12,20};
  int centmax[] = {4,12,20,40};
  
  float ajmin[] = { 0.00, 0.13, 0.24, 0.35 };
  float ajmax[] = { 0.13, 0.24, 0.35, 1.00 };
 
  
  TFile * outf = new TFile(Form("stdmcv2_%d.root",filenum),"recreate");
  
  for(int i = 0 ; i < 3 ; ++i)
  {
    cout<<"pt iteration "<<i<<endl;
    for(int cent = 0 ; cent < 4 ; ++cent)
    {
      cout<<"cent iteration "<<cent<<endl;
      for(int aj = 0 ; aj < 4 ; ++aj)
      {
        cout<<"aj iteration "<<aj<<endl;
        
        TH2D * ljtsig = JetTrackSignal    (0, leadingjetptlow[i], leadingjetpthigh[i] , subleadingjetptlow[i] , subleadingjetpthigh[i] , ptasslow[i] , ptasshigh[i], centmin[cent], centmax[cent],ajmin[aj],ajmax[aj]);
        TH2D * ljtbak = JetTrackBackground(0, leadingjetptlow[i], leadingjetpthigh[i] , subleadingjetptlow[i] , subleadingjetpthigh[i] , ptasslow[i] , ptasshigh[i], centmin[cent], centmax[cent],ajmin[aj],ajmax[aj]);
        TH2D * ljtcorr = (TH2D*)ljtsig->Clone(Form("corr_leadingjet%d_%d_ass%d_%d_cmin%d_cmax%d_ajmin%2.2f_ajmax%2.2f",(int)leadingjetptlow[i],(int)leadingjetpthigh[i],(int)ptasslow[i],(int)ptasshigh[i],centmin[cent],centmax[cent],ajmin[aj],ajmax[aj]));
        ljtcorr->Divide(ljtbak);
        ljtcorr->Scale(ljtbak->GetBinContent(ljtbak->FindBin(0,0)));
        ljtcorr->GetXaxis()->SetRange(ljtcorr->GetXaxis()->FindBin(-1.6),ljtcorr->GetXaxis()->FindBin(1.6));
        ljtcorr->GetYaxis()->SetRange(ljtcorr->GetYaxis()->FindBin(-3.1415926/2.0),ljtcorr->GetYaxis()->FindBin(3*3.1415926/2.0));
        
        TH2D * sljtsig = JetTrackSignal    (1, leadingjetptlow[i], leadingjetpthigh[i] , subleadingjetptlow[i] , subleadingjetpthigh[i] , ptasslow[i] , ptasshigh[i], centmin[cent], centmax[cent],ajmin[aj],ajmax[aj]);
        TH2D * sljtbak = JetTrackBackground(1, leadingjetptlow[i], leadingjetpthigh[i] , subleadingjetptlow[i] , subleadingjetpthigh[i] , ptasslow[i] , ptasshigh[i], centmin[cent], centmax[cent],ajmin[aj],ajmax[aj]);
        TH2D * sljtcorr = (TH2D*)sljtsig->Clone(Form("corr_subleadingjet%d_%d_ass%d_%d_cmin%d_cmax%d_ajmin%2.2f_ajmax%2.2f",(int)leadingjetptlow[i],(int)leadingjetpthigh[i],(int)ptasslow[i],(int)ptasshigh[i],centmin[cent],centmax[cent],ajmin[aj],ajmax[aj]));
        sljtcorr->Divide(sljtbak);
        sljtcorr->Scale(sljtbak->GetBinContent(sljtbak->FindBin(0,0)));
        sljtcorr->GetXaxis()->SetRange(sljtcorr->GetXaxis()->FindBin(-4.0),sljtcorr->GetXaxis()->FindBin(4.0));
        sljtcorr->GetYaxis()->SetRange(sljtcorr->GetYaxis()->FindBin(-3.1415926/2.0),sljtcorr->GetYaxis()->FindBin(3*3.1415926/2.0));
      }
    }
  }
  outf->Write();
  outf->Close();
}