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

void frameViewer(char* arg) {
	//Take the arguments and save them into respective strings
	std::string infileName;
	std::string inF;
	std::string inPrefix;
	std::string runs, layers;
	std::string runCount;
	std::istringstream stm(arg);

	inPrefix = "/home/p180f/Do_Jo_Ol_Ma/Analysis/MainProcedure/testMain/rawRoot/";

	while (true) {	
		if (std::getline(stm, layers, ' ')) {

			//load the input root files
			TChain *chain = new TChain("fourChamTree");

			for (int i=0; ; i++) {
				runCount = std::to_string(i);
				inF = "run" + runCount + "_" + layers + "layers.root";
				infileName = inPrefix + inF;
				ifstream fin;
				fin.open(infileName.c_str());
				if (!fin.fail()) {
					fin.close();
					chain->Add(infileName.c_str());
					std::cout << "Got " << inF << std::endl;
				} else break;
			}

			const int width=480;	//width of the raw image
			const int height=640;	//height of the raw image


			int x=-10;	//x from file
			int y=-10;	//y from file
			int intensity=-10;	//pixle intensity from file
			int pNum=0;//a counter of the order in which the frame was processed

			//the 2d array which will store each frame of image data.
			int frame[480][640]={0};

			//variables
			int UNIXtime=0;
			float tdc[2]={-10,-10};


			//TTree *T = new TTree("T","TTree of muplus data");
			//add the 'branches' to the tree we will now read in
			chain->SetBranchAddress("pNum",&pNum);	//branch for the frame number
			chain->SetBranchAddress("frame",&frame);	//branch for frame data

			TH2I *frameHisto = new TH2I("Single 4 Spark Event","",width,0,width,height,0,height);	//histogram for the stacked images

			TH1I *chamber1 = new TH1I("chamber1","Chamber 1",width,0,width);//histogram for chamber 1 (the top one)
			TH1I *chamber2 = new TH1I("chamber2","Chamber 2",width,0,width);//histogram for chamber 2
			TH1I *chamber3 = new TH1I("chamber3","Chamber 3",width,0,width);//histogram for chamber 3
			TH1I *chamber4 = new TH1I("chamber4","Chamber 4",width,0,width);//histogram for chamber 4 (the bottom one)
			
			TH1I *chamber1y = new TH1I("chamber1y","Chamber 1",height,0,height);//histogram for chamber 1 (the top one)
			TH1I *chamber2y = new TH1I("chamber2y","Chamber 2",height,0,height);//histogram for chamber 2
			TH1I *chamber3y = new TH1I("chamber3y","Chamber 3",height,0,height);//histogram for chamber 3
			TH1I *chamber4y = new TH1I("chamber4y","Chamber 4",height,0,height);//histogram for chamber 4 (the bottom one)

			//output the plot of the stacked images
		
			TCanvas *fH2 = new TCanvas("fH2", "Single 4 Spark Event", 0, 0, 800, 800);
			fH2->cd();
			frameHisto->Draw();
			frameHisto->GetXaxis()->SetTitle("X position (px)");
			frameHisto->GetXaxis()->CenterTitle();

			frameHisto->GetYaxis()->SetTitle("Intensity");
			frameHisto->GetYaxis()->SetTitleOffset(1.4);
			frameHisto->GetYaxis()->CenterTitle();

			TCanvas *pc2 = new TCanvas("pc2","Frame",0,0,800,800);
			
			pc2->Divide(2,2);
			
			pc2->cd(1);
			chamber1->Draw();
			chamber1->GetXaxis()->SetTitle("X position (px)");
			chamber1->GetXaxis()->CenterTitle();

			chamber1->GetYaxis()->SetTitle("Intensity");
			chamber1->GetYaxis()->SetTitleOffset(1.4);
			chamber1->GetYaxis()->CenterTitle();

			pc2->cd(2);
			chamber2->Draw();
			chamber2->GetXaxis()->SetTitle("X position (px)");
			chamber2->GetXaxis()->CenterTitle();

			chamber2->GetYaxis()->SetTitle("Intensity");
			chamber2->GetYaxis()->SetTitleOffset(1.4);
			chamber2->GetYaxis()->CenterTitle();

			pc2->cd(3);
			chamber3->Draw();
//.........这里部分代码省略.........

void CsIProj()
{
  TFile *file = new TFile("../root/NZ_55_New.root");
  TFile *gates = new TFile("../gates/zlines.root");
  TFile *gates2 = new TFile("../gates/zlines_new.root");

  ofstream ofile("CsI_55A_New.dat");
  
  TCanvas *mycan = (TCanvas*)gROOT->FindObjectAny("mycan");
  if(!mycan)
    {
      mycan = new TCanvas("mycan","mycan");
      mycan->Divide(1,2);
    }
  
  
  ostringstream outstring;
  string name;
  int p1= 30, p2=50; //+- fit limits up to 2 peaks. May be different.
  int const num_par = 5; //number of peaks times 2(pol1)+3(gaus).
  
  for(int ic =0;ic<56;ic++)
    {
      
      outstring.str("");
      outstring << "dEE/dEE_" << ic;
      name = outstring.str();
      
      mycan->cd(1);
      TH2I *hist = (TH2I*)file->Get(name.c_str());
      hist->Draw("col");  
      hist->GetXaxis()->SetRangeUser(200.,1800.); 
      hist->GetYaxis()->SetRangeUser(5.,50.); 
   
      
      if(ic <16 || ic > 31)
	TCutG *mycut = (TCutG*)gates->Get(Form("Zline_%i_2_4",ic));
      else 
	TCutG *mycut = (TCutG*)gates2->Get(Form("Zline_%i_2_4",ic));
      mycut->Draw();

      file->cd();
      outstring.str("");
      outstring << "CsI/CsIGate/ECsI_" << ic << "_Gate";
      name = outstring.str();
      gPad->SetLogz();

      mycan->cd(2);
      TH1I * proj = (TH1I*)file->Get(name.c_str());
      proj->Draw();
      proj->Rebin(4);
      proj->GetXaxis()->SetRangeUser(700.,1800.);

      mycan->Modified();
      mycan->Update();

      TMarker * mark;
      mark=(TMarker*)mycan->WaitPrimitive("TMarker"); //Get the Background limits
      int bkg_lo = mark->GetX();
      delete mark;  
      mark=(TMarker*)mycan->WaitPrimitive("TMarker");
      int bkg_hi = mark->GetX();
      delete mark;
      mark=(TMarker*)mycan->WaitPrimitive("TMarker"); // Get the 1st peak initial guess
      int peak1 = mark->GetX();
      delete mark;
      
      
      double par[num_par] = {0.};
      double out[num_par] = {0.}; 
      int peak1_lo = peak1 - p1, peak1_hi = peak1 + p1; // Peak center and limits
      
      
      TF1 *l1 = new TF1("l1", "pol1", bkg_lo, bkg_hi);
      TF1 *g1 = new TF1("g1", "gaus", peak1_lo,peak1_hi);
      
      TF1 *total = new TF1("total", "pol1(0)+gaus(2)", bkg_lo,bkg_hi);
      
      proj->Fit(l1,"R");
      proj->Fit(g1,"R+");
      
      l1->GetParameters(&par[0]);
      g1->GetParameters(&par[2]);
      
      total->SetParameters(par);
      proj->Fit(total,"R");
      total->GetParameters(out);
      
      
      ofile << ic << " " << out[3] << endl;
      
      outstring.str("");
      outstring << "55A_" << ic;
      name = outstring.str();
      total->SetName(name.c_str());
      total->Draw("same");
      mycan->Modified();
      mycan->Update();
      
      bool IsGood = 0;
//.........这里部分代码省略.........

int main (int argc, char** argv)
{

  if (argc < 3) {
    printHelp() ;
    exit (1) ;
  }

  std::string inputfiles, inputdir ;
  std::string outputRootName = "histoTPG.root" ;
  int verbose = 0 ;
  int occupancyCut = 0 ;
  std::string l1algo ; 

  bool ok(false) ;
  for (int i=0 ; i<argc ; i++) {
    if (argv[i] == std::string("-h") ) {
      printHelp() ;
      exit(1);
    }
    if (argv[i] == std::string("-i") && argc>i+1) {
      ok = true ;
      inputfiles = argv[i+1] ;
    }
    if (argv[i] == std::string("-d") && argc>i+1) inputdir = argv[i+1] ;
    if (argv[i] == std::string("-o") && argc>i+1) outputRootName = argv[i+1] ;
    if (argv[i] == std::string("-v") && argc>i+1) verbose = atoi(argv[i+1]) ;
    if (argv[i] == std::string("-l1") && argc>i+1) l1algo =  std::string(argv[i+1]) ;
    if (argv[i] == std::string("--cutTPOccup") && argc>i+1) occupancyCut = atoi(argv[i+1]) ;
  }
  if (!ok) {
    std::cout<<"No input files have been given: nothing to do!"<<std::endl ;
    printHelp() ;
    exit(1);
  }
  
  std::vector<int> algobits ;
  std::vector<std::string> algos = split(l1algo,",") ;
  for (unsigned int i=0 ; i<algos.size() ; i++) algobits.push_back(atoi(algos[i].c_str())) ;


  unsigned int ref = 2 ;



  ///////////////////////
  // book the histograms
  ///////////////////////

  TH2F * occupancyTP = new TH2F("occupancyTP", "Occupancy TP data", 72, 1, 73, 38, -19, 19) ;
  occupancyTP->GetYaxis()->SetTitle("eta index") ;
  occupancyTP->GetXaxis()->SetTitle("phi index") ;
  TH2F * occupancyTPEmul = new TH2F("occupancyTPEmul", "Occupancy TP emulator", 72, 1, 73, 38, -19, 19) ;
  occupancyTPEmul->GetYaxis()->SetTitle("eta index") ;
  occupancyTPEmul->GetXaxis()->SetTitle("phi index") ;

  TH1F * TP = new TH1F("TP", "TP", 256, 0., 256.) ;
  TP->GetXaxis()->SetTitle("TP (ADC)") ;
  TH1F * TPEmul = new TH1F("TPEmul", "TP Emulator", 256, 0., 256.) ;
  TPEmul->GetXaxis()->SetTitle("TP (ADC)") ;
  TH1F * TPEmulMax = new TH1F("TPEmulMax", "TP Emulator max", 256, 0., 256.) ;
  TPEmulMax->GetXaxis()->SetTitle("TP (ADC)") ;
  TH3F * TPspectrumMap3D = new TH3F("TPspectrumMap3D", "TP data spectrum map", 72, 1, 73, 38, -19, 19, 256, 0., 256.) ;
  TPspectrumMap3D->GetYaxis()->SetTitle("eta index") ;
  TPspectrumMap3D->GetXaxis()->SetTitle("phi index") ;

  TH1F * TPMatchEmul = new TH1F("TPMatchEmul", "TP data matching Emulator", 7, -1., 6.) ;
  TH1F * TPEmulMaxIndex = new TH1F("TPEmulMaxIndex", "Index of the max TP from Emulator", 7, -1., 6.) ;
  TH3I * TPMatchEmul3D = new TH3I("TPMatchEmul3D", "TP data matching Emulator", 72, 1, 73, 38, -19, 19, 7, -1, 6) ;
  TPMatchEmul3D->GetYaxis()->SetTitle("eta index") ;
  TPMatchEmul3D->GetXaxis()->SetTitle("phi index") ;

  TH2I * ttfMismatch = new TH2I("ttfMismatch", "TTF mismatch map",  72, 1, 73, 38, -19, 19) ;
  ttfMismatch->GetYaxis()->SetTitle("eta index") ;
  ttfMismatch->GetXaxis()->SetTitle("phi index") ;

  ///////////////////////
  // Chain the trees:
  ///////////////////////

  TChain * chain = new TChain ("EcalTPGAnalysis") ;
  std::vector<std::string> files ;
  if (inputfiles.find(std::string(",")) != std::string::npos) files = split(inputfiles,",") ;
  if (inputfiles.find(std::string(":")) != std::string::npos) {
    std::vector<std::string> filesbase = split(inputfiles,":") ;
    if (filesbase.size() == 4) {
      int first = atoi(filesbase[1].c_str()) ;
      int last = atoi(filesbase[2].c_str()) ;
      for (int i=first ; i<=last ; i++) {
	std::stringstream name ;
	name<<filesbase[0]<<i<<filesbase[3] ;
	files.push_back(name.str()) ;
      }
    }
  }
  for (unsigned int i=0 ; i<files.size() ; i++) {
    files[i] = inputdir+"/"+files[i] ;
    std::cout<<"Input file: "<<files[i]<<std::endl ;
    chain->Add (files[i].c_str()) ;
  }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
						if(frame[x][y]>0){
							frameHisto->Fill(x,y,frame[x][y]);
							frameHistoAll->Fill(x,y,frame[x][y]);
							if(y>580 && y<610){
								chamber1->Fill(x,frame[x][y]);
								chamber1All->Fill(x,frame[x][y]);
							}
							else if(y>400 && y<440){
								chamber2->Fill(x,frame[x][y]);
								chamber2All->Fill(x,frame[x][y]);
							}
							else if(y>240 && y<280){
								chamber3->Fill(x,frame[x][y]);
								chamber3All->Fill(x,frame[x][y]);
							}
							else if(y>50 && y<100){
								chamber4->Fill(x,frame[x][y]);
								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);

void plot_pad_size_in_layer(TString digiPar="trd.v13/trd_v13g.digi.par", Int_t nlines=1, Int_t nrows_in_sec=0, Int_t alllayers=1) 
{

  gStyle->SetPalette(1,0);
  gROOT->SetStyle("Plain");
  gStyle->SetPadTickX(1);                        
  gStyle->SetPadTickY(1); 
  gStyle->SetOptStat(kFALSE);
  gStyle->SetOptTitle(kFALSE);

  Bool_t read = false;
  TH2I *fLayerDummy = new TH2I("LayerDummy","",1200,-600,600,1000,-500,500);
  fLayerDummy->SetXTitle("x-coordinate [cm]");
  fLayerDummy->SetYTitle("y-coordinate [cm]");
  fLayerDummy->GetXaxis()->SetLabelSize(0.02);
  fLayerDummy->GetYaxis()->SetLabelSize(0.02);
  fLayerDummy->GetZaxis()->SetLabelSize(0.02);
  fLayerDummy->GetXaxis()->SetTitleSize(0.02);
  fLayerDummy->GetXaxis()->SetTitleOffset(1.5);
  fLayerDummy->GetYaxis()->SetTitleSize(0.02);
  fLayerDummy->GetYaxis()->SetTitleOffset(2);
  fLayerDummy->GetZaxis()->SetTitleSize(0.02);
  fLayerDummy->GetZaxis()->SetTitleOffset(-2);

  TString title;
  TString title1, title2, title3;
  TString buffer;
  TString firstModule = "";
  Int_t blockCounter(0), startCounter(0); // , stopCounter(0);
  Double_t msX(0), msY(0), mpX(0), mpY(0), mpZ(0), psX(0), psY(0);
  Double_t ps1X(0), ps1Y(0), ps2X(0), ps2Y(0), ps3X(0), ps3Y(0);
  Int_t modId(0), layerId(0);
  Double_t sec1(0), sec2(0), sec3(0);
  Double_t row1(0), row2(0), row3(0);
  std::map<float, TCanvas*> layerView;// map key is z-position of modules
  std::map<float, TCanvas*>::iterator it;
  ifstream digipar;

  digipar.open(digiPar.Data(), ifstream::in);
  while (digipar.good()) {
    digipar >> buffer;
    //cout << "(" << blockCounter << ")    " << buffer << endl;
    if (blockCounter == 19)
      firstModule = buffer;
    if (buffer == (firstModule + ":")){
      //cout << buffer << " <===========================================" << endl;
      read = true;
    }
    if (read) {
      startCounter++;
      if (startCounter == 1)   // position of module position in x
      {
	modId = buffer.Atoi();
        layerId = (modId & (15 << 4)) >> 4;  // from CbmTrdAddress.h
      }

      if (startCounter == 5)   // position of module position in x
	mpX = buffer.Atof();
      if (startCounter == 6)   // position of module position in y
	mpY = buffer.Atof();
      if (startCounter == 7)   // position of module position in z
	mpZ = buffer.Atof();
      if (startCounter == 8)   // position of module size in x
	msX = buffer.Atof();
      if (startCounter == 9)   // position of module size in y
	msY = buffer.Atof();

      if (startCounter == 12)   // sector 1 size in y
	sec1 = buffer.Atof();
      if (startCounter == 13)   // position of pad size in x - do not take the backslash (@14)
	ps1X = buffer.Atof();
      if (startCounter == 15)   // position of pad size in y
	ps1Y = buffer.Atof();

      if (startCounter == 17)   // sector 2 size in y
	sec2 = buffer.Atof();
      if (startCounter == 18)   // position of pad size in x
      {
	ps2X = buffer.Atof();
        psX = ps2X;   // for backwards compatibility - sector 2 is default sector
      }
      if (startCounter == 19)   // position of pad size in y
      {
	ps2Y = buffer.Atof();
        psY = ps2Y;   // for backwards compatibility - sector 2 is default sector
      }

      if (startCounter == 21)   // sector 3 size in y
	sec3 = buffer.Atof();
      if (startCounter == 22)   // position of pad size in x
	ps3X = buffer.Atof();
      if (startCounter == 23)   // position of pad size in y
	ps3Y = buffer.Atof();

//      if (startCounter == 23)   // last element
//      {
//        printf("moduleId         : %d, %d\n", modId, layerId);
//        printf("pad size sector 1: (%.2f cm, %.2f cm) pad area: %.2f cm2\n", ps1X, ps1Y, ps1X*ps1Y);
//        printf("pad size sector 2: (%.2f cm, %.2f cm) pad area: %.2f cm2\n", ps2X, ps2Y, ps2X*ps2Y);
//        printf("pad size sector 3: (%.2f cm, %.2f cm) pad area: %.2f cm2\n", ps3X, ps3Y, ps3X*ps3Y);
//.........这里部分代码省略.........

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