C++ TH1I::GetBinCenter方法代码示例

void ExtractTrackBasedTiming(TString fileName = "hd_root.root", int runNumber = 10390, TString variation = "default", bool verbose = false,TString prefix = ""){

   // set "prefix" in case you want to ship the txt files elsewhere...
   cout << "Performing Track Matched timing fits for File: " << fileName.Data() << " Run: " << runNumber << " Variation: " << variation.Data() << endl;

   ExtractTrackBasedTimingNS::thisFile = TFile::Open( fileName , "UPDATE");
   if (ExtractTrackBasedTimingNS::thisFile == 0) {
      cout << "Unable to open file " << fileName.Data() << "...Exiting" << endl;
      return;
   }

   //We need the existing constants, The best we can do here is just read them from the file.
   vector<double> sc_tdc_time_offsets;
   vector<double> sc_fadc_time_offsets;
   vector<double> tof_tdc_time_offsets;
   vector<double> tof_fadc_time_offsets;
   vector<double> tagm_tdc_time_offsets;
   vector<double> tagm_fadc_time_offsets;
   vector<double> tagh_tdc_time_offsets;
   vector<double> tagh_fadc_time_offsets;
   vector<double> tagh_counter_quality;

   double sc_t_base_fadc, sc_t_base_tdc;
   double tof_t_base_fadc, tof_t_base_tdc;
   double bcal_t_base_fadc, bcal_t_base_tdc;
   double tagm_t_base_fadc, tagm_t_base_tdc;
   double tagh_t_base_fadc, tagh_t_base_tdc;
   double fdc_t_base_fadc, fdc_t_base_tdc;
   double fcal_t_base;
   double cdc_t_base;
   double RF_Period;

   cout << "Grabbing CCDB constants..." << endl;
   // Base times
   GetCCDBConstants1("/CDC/base_time_offset" ,runNumber, variation, cdc_t_base);
   GetCCDBConstants1("/FCAL/base_time_offset",runNumber, variation, fcal_t_base);
   GetCCDBConstants1("/PHOTON_BEAM/RF/beam_period",runNumber, variation, RF_Period);
   GetCCDBConstants2("/FDC/base_time_offset" ,runNumber, variation, fdc_t_base_fadc, fdc_t_base_tdc);
   GetCCDBConstants2("/BCAL/base_time_offset" ,runNumber, variation, bcal_t_base_fadc, bcal_t_base_tdc);
   GetCCDBConstants2("/PHOTON_BEAM/microscope/base_time_offset" ,runNumber, variation, tagm_t_base_fadc, tagm_t_base_tdc);
   GetCCDBConstants2("/PHOTON_BEAM/hodoscope/base_time_offset" ,runNumber, variation, tagh_t_base_fadc, tagh_t_base_tdc);
   GetCCDBConstants2("/START_COUNTER/base_time_offset" ,runNumber, variation, sc_t_base_fadc, sc_t_base_tdc);
   GetCCDBConstants2("/TOF/base_time_offset" ,runNumber, variation, tof_t_base_fadc, tof_t_base_tdc);
   // Per channel
   //GetCCDBConstants("/BCAL/TDC_offsets"    ,runNumber, variation, bcal_tdc_offsets);
   //GetCCDBConstants("/FCAL/timing_offsets" ,runNumber, variation, fcal_adc_offsets);
   GetCCDBConstants("/START_COUNTER/adc_timing_offsets" ,runNumber, variation, sc_fadc_time_offsets);
   GetCCDBConstants("/START_COUNTER/tdc_timing_offsets" ,runNumber, variation, sc_tdc_time_offsets);
   GetCCDBConstants("/PHOTON_BEAM/microscope/fadc_time_offsets" ,runNumber, variation, tagm_fadc_time_offsets,3);// Interested in 3rd column
   GetCCDBConstants("/PHOTON_BEAM/microscope/tdc_time_offsets"  ,runNumber, variation, tagm_tdc_time_offsets,3);
   GetCCDBConstants("/PHOTON_BEAM/hodoscope/fadc_time_offsets"  ,runNumber, variation, tagh_fadc_time_offsets,2);// Interested in 2nd column
   GetCCDBConstants("/PHOTON_BEAM/hodoscope/tdc_time_offsets"   ,runNumber, variation, tagh_tdc_time_offsets,2);
   GetCCDBConstants("/PHOTON_BEAM/hodoscope/counter_quality"    ,runNumber, variation, tagh_counter_quality,2);
   GetCCDBConstants("/TOF/adc_timing_offsets",runNumber, variation, tof_fadc_time_offsets);
   GetCCDBConstants("/TOF/timing_offsets",runNumber, variation, tof_tdc_time_offsets);

   cout << "CDC base times = " << cdc_t_base << endl;
   cout << "FCAL base times = " << fcal_t_base << endl;
   cout << "FDC base times = " << fdc_t_base_fadc << ", " << fdc_t_base_tdc << endl;
   cout << "BCAL base times = " << bcal_t_base_fadc << ", " << bcal_t_base_tdc << endl;
   cout << "SC base times = " << sc_t_base_fadc << ", " << sc_t_base_tdc << endl;
   cout << "TOF base times = " << tof_t_base_fadc << ", " << tof_t_base_tdc << endl;
   cout << "TAGH base times = " << tagh_t_base_fadc << ", " << tagh_t_base_tdc << endl;
   cout << "TAGM base times = " << tagm_t_base_fadc << ", " << tagm_t_base_tdc << endl;

   cout << endl;
   cout << "RF_Period = " << RF_Period << endl;
   cout << endl;

   cout << "Done grabbing CCDB constants...Entering fits..." << endl;

   // Do our final step in the timing alignment with tracking

   //When the RF is present we can try to simply pick out the correct beam bucket for each of the runs
   //First just a simple check to see if we have the appropriate data
   bool useRF = false;
   TH1I *testHist = ExtractTrackBasedTimingNS::Get1DHistogram("HLDetectorTiming", "TAGH_TDC_RF_Compare","Counter ID 001");
   if (testHist != NULL){ // Not great since we rely on channel 1 working, but can be craftier later.
      cout << "Using RF Times for Calibration" << endl;
      useRF = true;
   }
   ofstream outFile;
   TH2I *thisHist; 
   thisHist = ExtractTrackBasedTimingNS::Get2DHistogram("HLDetectorTiming", "TRACKING", "TAGM - SC Target Time");
   if (useRF) thisHist = ExtractTrackBasedTimingNS::Get2DHistogram("HLDetectorTiming", "TRACKING", "TAGM - RFBunch Time");
   if (thisHist != NULL){
      //Statistics on these histograms are really quite low we will have to rebin and do some interpolation
      outFile.open(prefix + "tagm_tdc_timing_offsets.txt", ios::out | ios::trunc);
      outFile.close(); // clear file
      outFile.open(prefix + "tagm_adc_timing_offsets.txt", ios::out | ios::trunc);
      outFile.close(); // clear file
      int nBinsX = thisHist->GetNbinsX();
      int nBinsY = thisHist->GetNbinsY();
      TH1D * selectedTAGMOffset = new TH1D("selectedTAGMOffset", "Selected TAGM Offset; Column; Offset [ns]", nBinsX, 0.5, nBinsX + 0.5);
      TH1I * TAGMOffsetDistribution = new TH1I("TAGMOffsetDistribution", "TAGM Offset; TAGM Offset [ns]; Entries", 500, -250, 250);
      for (int i = 1 ; i <= nBinsX; i++){ 
         TH1D *projY = thisHist->ProjectionY("temp", i, i);
         // Scan over the histogram
         //chose the correct number of bins based on the histogram
         float nsPerBin = (projY->GetBinCenter(projY->GetNbinsX()) - projY->GetBinCenter(1)) / projY->GetNbinsX();
//.........这里部分代码省略.........

int main()
{

  float const conS =.3/2.35;
  ofstream fout("cal/backN.cal");
  ofstream fwhm("cal/fwhmback.dat");
  TFile f("sort.root");
  TCanvas* canvas[14];
  int Ntele = 14;
  int Nstrip = 32;
  ostringstream outstring;
  TH2I frame("frame","",10,4.5,9.5,10,0,130);
  frame.SetStats(kFALSE);

  double xx[14*32];
  double yy[14*32];

  TF1 *func = new TF1("fit",ThPeaks,3,9,4);
  double para[5];
  ifstream file("cal/back.cal");
  float intercept, slope;
  int i1,i2;
  string name;

  TH1F con("con","",500,0,10);
  for (int itele=0;itele<Ntele;itele++)
    {
      outstring.str("");
      outstring << "B"<<itele;
      name = outstring.str();
      canvas[itele] = new TCanvas(name.c_str());
      canvas[itele]->Divide(6,6);
      for (int istrip =0;istrip<Nstrip;istrip++)
        {
      
          canvas[itele]->cd(istrip+1);
          file >> i1 >> i2  >> slope >> intercept;


          outstring.str("");
          outstring << "back/cal/EBC"<<itele<<"_"<<istrip;
          string name = outstring.str();
          cout <<  name << endl;
          TH1I * hist = (TH1I*) f.Get(name.c_str());
	  frame.Draw();


          hist->SetStats(kFALSE);
          hist->GetXaxis()->SetRangeUser(4.5,9.5);
          con.GetXaxis()->SetRangeUser(4.5,9.5);
	  for (int i=1;i<=500;i++)
	    for (int j=1;j<500;j++)
	    {
              float deltax = hist->GetBinCenter(i)-con.GetBinCenter(j);
	      if (fabs(deltax) > 10.*conS)continue;
              float fact = gauss(deltax,0.,conS);
	      float y = fact*hist->GetBinContent(i)*hist->GetBinWidth(i);
              con.SetBinContent(j,y+con.GetBinContent(j));
	    }



	  for (int i=1;i<=500;i++) 
	    {
	      hist->SetBinContent(i,con.GetBinContent(i));
	      con.SetBinContent(i,0.);
	    }

          hist->Draw("same");



          func->SetParameter(0,0);
          func->SetParameter(1,1.);
	  func->FixParameter(2,conS);
          //func->SetParameter(2,0.1);
          func->SetParameter(3,8.);
          func->SetLineColor(2);
          //func->Draw("same");



	  hist->Fit(func);
          func->GetParameters(para);
          cout << "chisq=" << func->GetChisquare() << endl;
           if (fabs(para[1]-1.) < .2) 
	     { 

              slope *= para[1];
              intercept = intercept*para[1] + para[0];
	     }
            fout << itele << " " << istrip << " " 
                 << slope << " " << intercept << endl;
            fwhm << itele << " " << istrip << " " 
                  << para[2]*2.35 << endl;
            int ii = itele*32+istrip;
            xx[ii] = (float)ii;
            yy[ii] = para[2]*2.35;
            cout << para[0] << " " << para[1] << " " << para[2] << endl;

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

//.........这里部分代码省略.........
        }
    }
    inFile.close();


    // Do our final step in the timing alignment with tracking

    //When the RF is present we can try to simply pick out the correct beam bucket for each of the runs
    //First just a simple check to see if we have the appropriate data
    bool useRF = false;
    double RF_Period = 4.0080161;
    TH1I *testHist = Get1DHistogram("HLDetectorTiming", "TAGH_TDC_RF_Compare","Counter ID 001");
    if (testHist != NULL){ // Not great since we rely on channel 1 working, but can be craftier later.
        useRF = true;
    }
    ofstream outFile;
    TH2I *thisHist; 
    thisHist = Get2DHistogram("HLDetectorTiming", "TRACKING", "TAGM - SC Target Time");
    if (useRF) thisHist = Get2DHistogram("HLDetectorTiming", "TRACKING", "TAGM - RFBunch Time");
    if (thisHist != NULL){
        //Statistics on these histograms are really quite low we will have to rebin and do some interpolation
        outFile.open(prefix + "tagm_tdc_timing_offsets.txt", ios::out | ios::trunc);
        outFile.close(); // clear file
        outFile.open(prefix + "tagm_adc_timing_offsets.txt", ios::out | ios::trunc);
        outFile.close(); // clear file
        int nBinsX = thisHist->GetNbinsX();
        int nBinsY = thisHist->GetNbinsY();
        TH1D * selectedTAGMOffset = new TH1D("selectedTAGMOffset", "Selected TAGM Offset; Column; Offset [ns]", nBinsX, 0.5, nBinsX + 0.5);
        TH1I * TAGMOffsetDistribution = new TH1I("TAGMOffsetDistribution", "TAGM Offset; TAGM Offset [ns]; Entries", 500, -250, 250);
        for (int i = 1 ; i <= nBinsX; i++){ 
            TH1D *projY = thisHist->ProjectionY("temp", i, i);
            // Scan over the histogram
            //chose the correct number of bins based on the histogram
            float nsPerBin = (projY->GetBinCenter(projY->GetNbinsX()) - projY->GetBinCenter(1)) / projY->GetNbinsX();
            float timeWindow = 3; //ns (Full Width)
            int binWindow = int(timeWindow / nsPerBin);
            double maxEntries = 0;
            double maxMean = 0;
            for (int j = 1 ; j <= projY->GetNbinsX();j++){
                int minBin = j;
                int maxBin = (j + binWindow) <= projY->GetNbinsX() ? (j + binWindow) : projY->GetNbinsX();
                double sum = 0, nEntries = 0;
                for (int bin = minBin; bin <= maxBin; bin++){
                    sum += projY->GetBinContent(bin) * projY->GetBinCenter(bin);
                    nEntries += projY->GetBinContent(bin);
                    if (bin == maxBin){
                        if (nEntries > maxEntries) {
                            maxMean = sum / nEntries;
                            maxEntries = nEntries;
                        }
                    } 
                }
            }
            //In the case there is RF, our job is to pick just the number of the correct beam bunch, so that's really all we need.
            if(useRF) {
                int beamBucket = int((maxMean / RF_Period) + 0.5); // +0.5 to handle rounding correctly
                selectedTAGMOffset->SetBinContent(i, beamBucket);
                TAGMOffsetDistribution->Fill(beamBucket);
            }
            else{
                selectedTAGMOffset->SetBinContent(i, maxMean);
                TAGMOffsetDistribution->Fill(maxMean);
            }
        }
        /*
        if (!useRF){