C++ TH2I::ProjectionY方法代码示例

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();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
								chamber4All->Fill(x,frame[x][y]);
							}

						}

					}
				}
				cout << "Stacking frame number " << pNum << "\r";//this overwrites the line every time
			}
			cout << endl;

			//output the plot of the stacked images
			TCanvas *pc2 = new TCanvas("pc2","Stacked Frames",0,0,600,800);
			pc2->cd();
			frameHisto->SetStats(false);
			frameHisto->Draw("colz");
			frameHisto->GetXaxis()->SetTitle("X position (px)");
			//frameHisto->GetXaxis()->SetTitleSize(0.055);
			//frameHisto->GetXaxis()->SetTitleOffset(1.0);
			//frameHisto->GetXaxis()->SetLabelSize(0.055);
			frameHisto->GetXaxis()->CenterTitle();

			frameHisto->GetYaxis()->SetTitle("Y position (px)");
			//frameHisto->GetYaxis()->SetTitleSize(0.055);
			//frameHisto->GetYaxis()->SetTitleOffset(0.9);
			//frameHisto->GetYaxis()->SetLabelSize(0.055);
			frameHisto->GetYaxis()->CenterTitle();
			gPad->Update();
//			pc2->Print("chamberStack.png");//output to a graphics file 

			//plot the projection onto the Y axis (so we can find our cuts in Y to select each chamber)
			TCanvas *projC = new TCanvas("projC","",0,0,800,600);
			projC->cd();
			TH1D *ydist = frameHisto->ProjectionY("ydist");
			ydist->Draw();
			ydist->GetYaxis()->SetTitle("Entries");
			ydist->GetYaxis()->CenterTitle();

			TCanvas *sliceX = new TCanvas("sliceX","",0,0,800,600);

			sliceX->Divide(2,2);
			
			sliceX->cd(1);
			chamber1->Draw();
			chamber1->GetXaxis()->SetTitle("X position (px)");
			chamber1->GetXaxis()->CenterTitle();
			chamber1->GetYaxis()->SetTitle("Y position (px)");
			chamber1->GetYaxis()->CenterTitle();
		//	chamber1->GetYaxis()->SetMaxDigits(2);
			
			sliceX->cd(2);
			chamber2->Draw();
			chamber2->GetXaxis()->SetTitle("X position (px)");
			chamber2->GetXaxis()->CenterTitle();
			chamber2->GetYaxis()->SetTitle("Y position (px)");
			chamber2->GetYaxis()->CenterTitle();
		//	chamber2->GetYaxis()->SetMaxDigits(2);
			
			sliceX->cd(3);
			chamber3->Draw();
			chamber3->GetXaxis()->SetTitle("X position (px)");
			chamber3->GetXaxis()->CenterTitle();
			chamber3->GetYaxis()->SetTitle("Y position (px)");
			chamber3->GetYaxis()->CenterTitle();
		//	chamber3->GetYaxis()->SetMaxDigits(2);
			

//.........这里部分代码省略.........
        while (getline (inFile, line)){
            istringstream iss(line);
            iss>>cdc_t_base; 
        }
    }
    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);
            }

// Do the extraction of the actual constants
void ExtractCDCDeformation(TString filename = "hd_root.root"){

    // Open our input and output file
    thisFile = TFile::Open(filename);
    TFile *outputFile = TFile::Open("CDCDeformation_Results.root", "RECREATE");

    // Check to make sure it is open
    if (thisFile == 0) {
        cout << "Unable to open file " << filename.Data() << "...Exiting" << endl;
        return;
    }

    // This stream will be for outputting the results in a format suitable for the CCDB
    // Will wait to open until needed
    ofstream textFile;
    textFile.open("CDC_Deformation.txt");

    // We want to display the direction of the shift as well as the magnitude in the "CDC view"
    // Let's make it happen
    int straw_offset[29] = {0,0,42,84,138,192,258,324,404,484,577,670,776,882,1005,1128,1263,1398,1544,1690,1848,2006,2176,2346,2528,2710,2907,3104,3313};
    int Nstraws[28] = {42, 42, 54, 54, 66, 66, 80, 80, 93, 93, 106, 106, 123, 123, 135, 135, 146, 146, 158, 158, 170, 170, 182, 182, 197, 197, 209, 209};
    double radius[28] = {10.72134, 12.08024, 13.7795, 15.14602, 18.71726, 20.2438, 22.01672, 23.50008, 25.15616, 26.61158, 28.33624, 29.77388, 31.3817, 32.75838, 34.43478, 35.81146, 38.28542, 39.7002, 41.31564, 42.73042, 44.34078, 45.75302, 47.36084, 48.77054, 50.37582, 51.76012, 53.36286, 54.74716};
    double phi[28] = {0, 0.074707844, 0.038166294, 0.096247609, 0.05966371, 0.012001551, 0.040721951, 0.001334527, 0.014963808, 0.048683644, 0.002092645, 0.031681749, 0.040719354, 0.015197341, 0.006786058, 0.030005892, 0.019704045, -0.001782064, -0.001306618, 0.018592421, 0.003686784, 0.022132975, 0.019600866, 0.002343723, 0.021301449, 0.005348855, 0.005997358, 0.021018761};

    TH2D * Amplitude_view[29];
    TH2D * Direction_view[29];
    TH2D * Vertical_view[29];
    TH2D * Horizontal_view[29];

    outputFile->mkdir("PerRing");
    outputFile->cd("PerRing");
    for(unsigned int iring=0; iring<28; iring++){
        double r_start = radius[iring] - 0.8;
        double r_end = radius[iring] + 0.8;
        double phi_start = phi[iring]; 
        double phi_end = phi_start + TMath::TwoPi();

        char hname[256];
        sprintf(hname, "Amplitude_view_ring[%d]", iring+1);
        Amplitude_view[iring+1] = new TH2D(hname, "", Nstraws[iring], phi_start, phi_end, 1, r_start, r_end);
        sprintf(hname, "Direction_view_ring[%d]", iring+1);
        Direction_view[iring+1] = new TH2D(hname, "", Nstraws[iring], phi_start, phi_end, 1, r_start, r_end);
        sprintf(hname, "Vertical_view_ring[%d]", iring+1);
        Vertical_view[iring+1] = new TH2D(hname, "", Nstraws[iring], phi_start, phi_end, 1, r_start, r_end);
        sprintf(hname, "Horizontal_view_ring[%d]", iring+1);
        Horizontal_view[iring+1] = new TH2D(hname, "", Nstraws[iring], phi_start, phi_end, 1, r_start, r_end);
    }

    //Fit function for 
    TF1 *f1 = new TF1("f1", "[0] + [1] * TMath::Cos(x + [2])", -3.14, 3.14);
    f1->SetParLimits(0, 0.5, 1.0);
    f1->SetParLimits(1, 0.0, 0.35);
    //f1->SetParLimits(2, -3.14, 3.14);
    f1->SetParameters(0.78, 0.0, 0.0);

    outputFile->cd();
    outputFile->mkdir("FitParameters");
    outputFile->cd("FitParameters");

    // Make some histograms to get the distributions of the fit parameters
    TH1I *h1_c0 = new TH1I("h1_c0", "Distribution of Constant", 100, 0.5, 1.0);
    TH1I *h1_c1 = new TH1I("h1_c1", "Distribution of Amplitude", 100, 0.0, 0.35);
    TH1I *h1_c2 = new TH1I("h1_c2", "Direction of Longest Drift Time", 100, -3.14, 3.14);
    TH1F *h1_c2_weighted = new TH1F("h1_c2_weighted", "Distribution of Direction weighted by amplitude", 100, -3.14, 3.14);
    TH2I *h2_c0_c1 = new TH2I("h2_c0_c1", "c_{1} Vs. c_{0}; c_{0}; c_{1}", 100, 0.5, 1.0, 100, 0, 0.35);
    TH2I *h2_c0_c2 = new TH2I("h2_c0_c2", "c_{2} Vs. c_{0}; c_{0}; c_{2}", 100, 0.5, 1.0, 100, -10, 10);
    TH2I *h2_c1_c2 = new TH2I("h2_c1_c2", "c_{2} Vs. c_{1}; c_{1}; c_{2}", 100, 0.0, 0.35, 100, -10, 10);

    outputFile->cd();
    outputFile->mkdir("Fits");
    outputFile->cd("Fits");

    // Now we want to loop through all available module/layer/sector and try to make a fit of each one
    int ring = 1, straw = 1;
    while (ring <= 28){
        cout << "Entering Fit " << endl;
        char folder[100];
        sprintf(folder, "Ring %.2i", ring);
        char strawname[100];
        sprintf(strawname,"Straw %.3i Predicted Drift Distance Vs phi_DOCA", straw);
        TH2I *thisStrawHistogram = Get2DHistogram("CDC_Cosmic_Per_Straw",folder,strawname);

        if (thisStrawHistogram != NULL) {

            // Now to do our fits. This time we know there are 16 bins.
            double percentile95[16], percentile97[16], percentile99[16]; // Location of 95, 97,and 99th percentile bins
            double binCenter[16];
            char name[100];
            sprintf(name,"Ring %.2i Straw %.3i", ring, straw);
            TH1D *extractedPoints = new TH1D(name, name, 16, -3.14, 3.14);
            for (int i = 1; i <= thisStrawHistogram->GetNbinsX() ; i++){
                TH1D *projY = thisStrawHistogram->ProjectionY(" ", i, i);
                binCenter[i-1] = thisStrawHistogram->GetXaxis()->GetBinCenter(i);
                int nbins = projY->GetNbinsX();
                //Get the total nubmer of entries
                int nEntries = projY->GetEntries();
                if (nEntries == 0) continue;
                double errorFraction = TMath::Sqrt(nEntries) / nEntries;
                double perc95 = 0.95*nEntries, perc97 = 0.97 * nEntries, perc99 = 0.99 * nEntries;
//.........这里部分代码省略.........