C++ TRandom3::SetSeed方法代码示例

RooDataHist * genHistFromModelPdf(const char * name, RooAbsPdf *model,  RooRealVar *var,    double ScaleLumi,  int range, int rebin, int seed ) {
  double genEvents =  model->expectedEvents(*var);
  TRandom3 *rndm = new TRandom3();
  rndm->SetSeed(seed);
  double nEvt = rndm->PoissonD( genEvents) ;
  int intEvt = ( (nEvt- (int)nEvt) >= 0.5) ? (int)nEvt +1 : int(nEvt);
  RooDataSet * data = model->generate(*var ,   intEvt   );
  cout<< " expected events for " << name << " = "<< genEvents << endl; 
  cout<< " data->numEntries() for name " << name << " == " << data->numEntries()<< endl;
  // cout<< " nEvt from PoissonD for" << name << " == " << nEvt<< endl;
  //cout<< " cast of nEvt  for" << name << " == " << intEvt<< endl; 
  RooAbsData *binned_data = data->binnedClone();
  TH1 * toy_hist = binned_data->createHistogram( name, *var, Binning(range/rebin )  );
  for(int i = 1; i <= toy_hist->GetNbinsX(); ++i) {
    toy_hist->SetBinError( i,  sqrt( toy_hist->GetBinContent(i)) );
    if(toy_hist->GetBinContent(i) == 0.00) {
      cout<< " WARNING: histo " << name << " has 0 enter in bin number " << i << endl;   
    }
    if(toy_hist->GetBinContent(i) < 0.1) {
      toy_hist->SetBinContent(i, 0.0);
      toy_hist->SetBinError(i, 0.0);
      cout<< " WARNING: setting value 0.0 to histo " << name << " for bin number " << i << endl;   
    }  
  }
  RooDataHist * toy_rooHist = new RooDataHist(name, name , RooArgList(*var), toy_hist );
  return toy_rooHist; 
}

int main(int argc, char** argv)
{
    file = TFile::Open("~/lhcb/rd/Kll/tuples/fromPatrick/Kmm_Q17_reduced.root");
    DecayTree = dynamic_cast<TTree*>(file->Get("finalTree_KMuMu"));
    cosTheta = new RooRealVar("cosThetaL", "cosThetaL", -1., 1.);
    sWeight  = new RooRealVar("sWeight",   "sWeight",   -5., 5.);
    data = new RooDataSet("data", "", DecayTree, RooArgList(*cosTheta,*sWeight), 0, "sWeight");

    fout.open("results.txt",std::ofstream::out);

    r.SetSeed(123456);

    Double_t a(0.), b(0.);
    for(Int_t j=-2; j<3; ++j) {
        for(Int_t k=0; k<10; ++k) {
            for(Int_t i=0; i<20; ++i) {
                a = k/-10.;
                b = j/10.;
                fit(i,a,b);
            }
        }
    }

    fout.close();

    return 0;
}

void plot_distortion(){
    
    TRandom3* rand = new TRandom3();

//    TF1 *func = new TF1("func","TMath::Abs([0]+[1]*x)",0,12);
//
//    rand->SetSeed(0);
//    Double_t slope=0.2*rand->Gaus(0,1);
//    Double_t anchor_point=3.5;
//    //want offset to be set by requiring func at anchor point to be 1
//    Double_t offset=(1-slope*anchor_point);
//    func->SetParameter(0,offset);
//    func->SetParameter(1,slope);
    
    //FN distortion
    TF1 *func = new TF1("func","[0]/(0.2*pow(x,0.1))+[1]",0,12);
    rand->SetSeed(0);
    Double_t scaling =0.2*rand->Gaus(0,1);
    func->SetParameter(0,scaling);
    func->SetParameter(1,0);
    //set offset so that func(FinalEnergy)=1;
    func->SetParameter(1,1-1*func->Eval(12));
    
    func->Draw();
    
    
   }

int decroissance_pi(float _p_pi = 1.0){

  //	random.seed(60934386)
  TRandom3 rand;
  rand.SetSeed();
  //float _p_pi;
  //cout << "Entrez l'impulsion des pions (en GeV) : p = ";
  //cin >> _p_pi;

  double x[1000] = {0.0};
  double y[1000] = {0.0};

  for (int i=0;i<1000;i++){
    //double _theta_cm_mu = rand.Uniform(TMath::Pi());
    double _theta_cm_mu = TMath::ACos(gRandom->Uniform(2.0)-1.0);
    double _phi_cm_mu = rand.Uniform(2*TMath::Pi());
    if (i < 10){
      cout << "i: " << i <<  " th_cm_mu = " << _theta_cm_mu << endl;
    }
    double _theta_lab_mu = 0.0, _p_lab_mu = 0.0;
    ThetaLab_mu(_p_pi, _theta_cm_mu, _theta_lab_mu, _p_lab_mu);
    cout << "th_lab_mu= " << _theta_lab_mu << " _p_lab_mu= " << _p_lab_mu << endl;
    x[i] = 1000*_theta_lab_mu;
    y[i] = _p_lab_mu;
  }

  TGraph *graph = new TGraph(1000, x, y);
  graph->SetTitle("Theta vs p_{LAB}");

  graph->Draw("A*");

  return 0;

}

void histo_graph(TH1D *histo,TGraph *graph,Int_t N_to_fill) {
  Int_t N_Points=graph->GetN();
  Double_t x,y;

  Double_t y_min=0.0;
  Double_t y_max=0.0;
  Double_t x_min=0.0;
  Double_t x_max=0.0;

  for (Int_t i=0;i<N_Points;i++) {
    graph->GetPoint(i,x,y);
    if (y>y_max) y_max=y;
    if (x>x_max) x_max=x;
    if (y<y_min) y_min=y;
    if (x<x_min) x_min=x;
  }
  Double_t histo_x_min=histo->GetXaxis()->GetXmin();
  Double_t histo_x_max=histo->GetXaxis()->GetXmax();
  if (histo_x_min>x_min) x_min=histo_x_min;
  if (histo_x_max<x_max) x_max=histo_x_max;

  //reverse graph order if necessary
  TGraph graph2(N_Points);
  Double_t x0,xNm1;
  graph->GetPoint(0,x0,y);
  graph->GetPoint(N_Points-1,xNm1,y);
  if (xNm1<x0) {
    for (Int_t iPoint=0;iPoint<N_Points;iPoint++) {
      Double_t x,y;
      graph->GetPoint(N_Points-iPoint,x,y);
      graph2.SetPoint(iPoint,x,y);
    }
    graph=&graph2;
  }

  TRandom3 r;
  r.SetSeed(0);

  Int_t N_tried=0;
  Int_t N_filled=0;
  while (N_filled<N_to_fill) {
    N_tried++;
    //rejection method
    x=r.Uniform(x_min,x_max);
    y=r.Uniform(0.0,y_max);
    
    if (y<graph->Eval(x)) {
      histo->Fill(x);
      N_filled++;
    }
  }
  printf("x_max=%8.6g; x_min=%8.6g; y_max=%8.6g; N_filled=%d; N_tried=%d\n",x_max,x_min,y_max,N_filled,N_tried);
  Double_t norm=(x_max-x_min)*y_max*N_filled/N_tried;
  histo->Scale(norm/histo->Integral("width"));

}

void createTBevents(int input){

  printf("Starting Simulation of data\n");

  //creating the output file
  char outputFileName[100] = {"OutputFile.root"};

  printf("Creating output file: %s \n",outputFileName);
  TFile * outputFile = new TFile(outputFileName,"RECREATE");


  //Counter for event number
  unsigned int eventNr;

  //Counter for total number of hits
  unsigned int hitsTotal = 0;

  
  short int col, row, adc;
  short int ladder = 2;
  short int mod = 3;
  short int disk = 2;
  short int blade = 2;
  short int panel = 2;
  
  //create the tree to store the data
  TTree *bpixTree[3];

  char title[30];
  for (int i=1; i<4; i++){
    sprintf(title,"BPIX_Digis_Layer%1d",i);
    bpixTree[i-1]= new TTree(title,title);
    bpixTree[i-1]->Branch("Event", &eventNr, "Event/i");           
    bpixTree[i-1]->Branch("Ladder", &ladder, "Ladder/S");      
    bpixTree[i-1]->Branch("Module", &mod, "Module/S");      
    bpixTree[i-1]->Branch("adc", &adc, "adc/S");       
    bpixTree[i-1]->Branch("col", &col, "col/S");       
    bpixTree[i-1]->Branch("row", &row, "row/S");
  }

  

  //Maximum number of events. Events does not correspond with Hits
  unsigned int maxEventNr = input;

  //Number of Hits per Event
  //This should be randomized later and be dependant on the rate
  double meanHitsPerEvent = 2;

  //Maximum particle flux [MHz cm^-2]
  int maxParticleFlux = 500;

  //number of hits in current event
  int hitsInEvent = -1; 

  //create a random number generator
  TRandom3 * random = new TRandom3();
  TRandom3 * randomrow = new TRandom3();
  TRandom3 * randomcol = new TRandom3();
  TRandom3 * randomadc = new TRandom3();

  //using custom function to distribute values
  //values used from http://ntucms1.cern.ch/pixel_dev/flux/v8/016031/fitspot_bin_11.pdf
  
  TF1 * fx = new TF1("xfunc", "[0]*exp(2.59349*exp(-0.5*((-x+3.24273/.15+[1]/.15)/7.07486*.15)**2)+2.07765*exp(-0.5*((-x+9.33060e-01/.15+[1]/.15)/2.24067*.15)**2)-4.21847)",0, 52);
  
  fx->SetParameters(1,337.0);
  fx->SetParameters(2,1.74);
  
  TF1 * fy = new TF1("yfunc", AsyGaus,0,80,4);
  
  fy->SetParNames("mean","sigma1","sigma2","amplitude");
  fy->SetParameter("mean",43);
  fy->SetParameter("sigma1",11.4);
  fy->SetParameter("sigma2",15.0);
  fy->SetParameter("amplitude",347.0);
  
  TF1 * fadc = new TF1("adcfunc", langaufun,0,400,4);
  fadc->SetParNames("scale","mpv","area","sigma");
  fadc->SetParameter("scale",19);
  fadc->SetParameter("mvp",220);
  fadc->SetParameter("area",10000);
  fadc->SetParameter("sigma",30);




  while (eventNr < maxEventNr){
    //printf("eventNr: %d \n",eventNr);
    //Function used for fitting according to Xin an Stefano
    
    //Start by generating the number of hits per event
    //following a poisson distribution
    random->SetSeed(0);
    hitsInEvent = random->Poisson(meanHitsPerEvent);
    //printf("hitsInEvent %d \n", hitsInEvent);

    hitsTotal += hitsInEvent;

    if(hitsInEvent < 0){
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  //   ue_fluctuation[3] = 0.1;
  //   ue_fluctuation[4] = 0.1;
  //   ue_fluctuation[5] = 0.01;
  // }
  // if(radius == 3) {
  //   ue_fluctuation[0] = 3.17;
  //   ue_fluctuation[1] = 0.7;
  //   ue_fluctuation[2] = 0.1;
  //   ue_fluctuation[3] = 0.1;
  //   ue_fluctuation[4] = 0.1;
  //   ue_fluctuation[5] = 0.01;
  // }
  // if(radius == 2) {
  //   ue_fluctuation[0] = 1.5;
  //   ue_fluctuation[1] = 0.4;
  //   ue_fluctuation[2] = 0.2;
  //   ue_fluctuation[3] = 0.1;
  //   ue_fluctuation[4] = 0.1;
  //   ue_fluctuation[5] = 0.01;
  // }

  // float  ue_fluctuation = 6.0;
  // if(radius == 4) ue_fluctuation = 6.0;
  // if(radius == 3) ue_fluctuation = 4.0;
  // if(radius == 2) ue_fluctuation = 2.0;

  // the UE Smear corresponds to doing recpt = recpt * (1 + rnd.Gaus(0, (float)ue_fluctuation[cBin]/recpt);
  
  // now start the event loop for each file. 
  if(printDebug) cout<<"Running through all the events now"<<endl;
  Long64_t nentries = jetpbpb[0]->GetEntries();
  if(printDebug) nentries = 10;
  TRandom3 rnd;
  rnd.SetSeed(endfile);

  for(int nEvt = 0; nEvt < nentries; ++ nEvt) {

    if(nEvt%10000 == 0)cout<<nEvt<<"/"<<nentries<<endl;
    if(printDebug)cout<<"nEvt = "<<nEvt<<endl;
    
    jetpbpb[0]->GetEntry(nEvt);
    jetpbpb[1]->GetEntry(nEvt);
    jetpbpb[2]->GetEntry(nEvt);
    jetpbpb[3]->GetEntry(nEvt);

    if(printDebug) cout<<"forest values = "<<hiBin_F<<", "<<evt_F<<", "<<run_F<<", "<<lumi_F<<", "<<vz_F<<endl;
    
    if(pHBHENoiseFilter_F == 0) continue; 
    if(pcollisionEventSelection_F == 0) continue; 
    if(fabs(vz_F)>15) continue;
    // if(!isGoodEvent_eS) continue; 

    int cBin = findBin(hiBin_F);//tells us the centrality of the event.
    if(cBin == -1) continue;

    // int jetCounter = 0;
    
    // for(int g = 0;g<nref_F;g++){
      
    //   if(eta_F[g]>=-2 && eta_F[g]<2){ //to select inside 
	
    // 	if(pt_F[g]>=50) jetCounter++;
	
    //   }//eta selection cut
      
    // }// jet loop

int main(int argc, char * argv[])
{
  if (argc!=3) {
    cerr << "Usage: " << argv[0] << " input_file output_file." << endl;
    return 1;
  }
  // variables in tree
  int run_id;
  int event_id;
  string * parent = 0;
  double t0;
  double t1;
  double energy;
  double maxDx;
  double maxDy;
  double maxDz;
  double maxDd;
  double primaryX;
  double primaryY;
  double primaryZ;

  // open the file for read
  TFile f(argv[1], "READ");
  TTree * tree = (TTree *) f.Get("simple_out");
  TBranch * b_parent;
  // mapping
  tree->SetMakeClass(1);
  tree->SetBranchAddress("runId", &run_id);
  tree->SetBranchAddress("eventId", &event_id);
  tree->SetBranchAddress("parent", &parent, &b_parent);
  tree->SetBranchAddress("t0", &t0);
  tree->SetBranchAddress("t1", &t1);
  tree->SetBranchAddress("energy", &energy);
  tree->SetBranchAddress("maxDx", &maxDx);
  tree->SetBranchAddress("maxDy", &maxDy);
  tree->SetBranchAddress("maxDz", &maxDz);
  tree->SetBranchAddress("maxDd", &maxDd);
  tree->SetBranchAddress("primaryX", &primaryX);
  tree->SetBranchAddress("primaryY", &primaryY);
  tree->SetBranchAddress("primaryZ", &primaryZ);

  // ROI
  double fwhm_0_5 = 0.005 * Q_value;
  double fwhm_1 = 0.01 * Q_value;
  double fwhm_3 = 0.03 * Q_value;

  double sigma_0_5 = getSigma(fwhm_0_5);
  double sigma_1 = getSigma(fwhm_1);
  double sigma_3 = getSigma(fwhm_3);

  double l_edge_0_5 = Q_value - 2 * sigma_0_5;
  double u_edge_0_5 = Q_value + 2 * sigma_0_5;
  double l_edge_1 = Q_value - 2 * sigma_1;
  double u_edge_1 = Q_value + 2 * sigma_1;
  double l_edge_3 = Q_value - 2 * sigma_3;
  double u_edge_3 = Q_value + 2 * sigma_3;
  cout << "0.5% FHWM " << fwhm_0_5 << " - sigma " << sigma_0_5 << endl;
  cout << "1.0% FHWM " << fwhm_1 << " - sigma " << sigma_1 << endl;
  cout << "3.0% FHWM " << fwhm_3 << " - sigma " << sigma_3 << endl;

  cout << "0.5% FHWM - (" << l_edge_0_5 << ", " << u_edge_0_5 << ")." << endl;
  cout << "1.0% FHWM - (" << l_edge_1 << ", " << u_edge_1 << ")." << endl;
  cout << "3.0% FHWM - (" << l_edge_3 << ", " << u_edge_3 << ")." << endl;

  // output file, only with the smeared energy.
  TFile fo(argv[2], "RECREATE");
  TTree * out_tree = new TTree ("smear_e", "smeared energy");
  // output variables
  double e_smear_0_5, e_smear_1, e_smear_3;
  string pparent;
  out_tree->Branch("runId", &run_id, "runId/I");
  out_tree->Branch("eventId", &event_id, "eventId/I");
  out_tree->Branch("parent", &pparent);
  out_tree->Branch("energy", &energy, "energy/D");
  out_tree->Branch("e_smear_0_5", &e_smear_0_5, "e_smear_0_5/D");
  out_tree->Branch("e_smear_1", &e_smear_1, "e_smear_1/D");
  out_tree->Branch("e_smear_3", &e_smear_3, "e_smear_3/D");
  out_tree->Branch("maxDx", &maxDx, "maxDx/D");
  out_tree->Branch("maxDy", &maxDy, "maxDy/D");
  out_tree->Branch("maxDz", &maxDz, "maxDz/D");
  out_tree->Branch("maxDd", &maxDd, "maxDd/D");
  out_tree->Branch("primaryX", &primaryX, "primaryX/D");
  out_tree->Branch("primaryY", &primaryY, "primaryY/D");
  out_tree->Branch("primaryZ", &primaryZ, "primaryZ/D");

  long nEntries = tree->GetEntries();
  cout << nEntries << endl;
  // loop over entries and smear the energy

  tr.SetSeed(0);
  // TCanvas c1("c1", "c1", 800, 600);
  // TH1D * th_esm = new TH1D("th_esm", "smeared energy", 500, 2200, 2700);
  // for (int i=0; i<1000; ++i) {
  //   double esm = smearEnergy(Q_value, sigma_3);
  //   th_esm->Fill(esm);
  // }
  // th_esm->Draw();
  // c1.Print("esm.png");
  vector<double> es;		// selected energy
  int n_0_5(0), ns_0_5(0);
//.........这里部分代码省略.........

void kdTreeBinning() {

   // -----------------------------------------------------------------------------------------------
   //  C r e a t e  r a n d o m  s a m p l e  w i t h  r e g u l a r  b i n n i n g  p l o t t i n g
   // -----------------------------------------------------------------------------------------------

   const UInt_t DATASZ = 10000;
   const UInt_t DATADIM = 2;
   const UInt_t NBINS = 50;

   Double_t smp[DATASZ * DATADIM];

   double mu[2] = {0,2};
   double sig[2] = {2,3};
   TRandom3 r;
   r.SetSeed(1);
   for (UInt_t i = 0; i < DATADIM; ++i)
      for (UInt_t j = 0; j < DATASZ; ++j)
         smp[DATASZ * i + j] = r.Gaus(mu[i], sig[i]);

   UInt_t h1bins = (UInt_t) sqrt(NBINS);

   TH2D* h1 = new TH2D("h1BinTest", "Regular binning", h1bins, -5., 5., h1bins, -5., 5.);
   for (UInt_t j = 0; j < DATASZ; ++j)
      h1->Fill(smp[j], smp[DATASZ + j]);


   // ---------------------------------------------------------------------------------------------
   // C r e a t e  K D T r e e B i n n i n g  o b j e c t  w i t h  T H 2 P o l y  p l o t t i n g
   // ---------------------------------------------------------------------------------------------

   TKDTreeBinning* kdBins = new TKDTreeBinning(DATASZ, DATADIM, smp, NBINS);

   UInt_t nbins = kdBins->GetNBins();
   UInt_t dim   = kdBins->GetDim();

   const Double_t* binsMinEdges = kdBins->GetBinsMinEdges();
   const Double_t* binsMaxEdges = kdBins->GetBinsMaxEdges();

   TH2Poly* h2pol = new TH2Poly("h2PolyBinTest", "KDTree binning", kdBins->GetDataMin(0), kdBins->GetDataMax(0), kdBins->GetDataMin(1), kdBins->GetDataMax(1));

   for (UInt_t i = 0; i < nbins; ++i) {
      UInt_t edgeDim = i * dim;
      h2pol->AddBin(binsMinEdges[edgeDim], binsMinEdges[edgeDim + 1], binsMaxEdges[edgeDim], binsMaxEdges[edgeDim + 1]);
   }

   for (UInt_t i = 1; i <= kdBins->GetNBins(); ++i)
      h2pol->SetBinContent(i, kdBins->GetBinDensity(i - 1));

   std::cout << "Bin with minimum density: " << kdBins->GetBinMinDensity() << std::endl;
   std::cout << "Bin with maximum density: " << kdBins->GetBinMaxDensity() << std::endl;

   TCanvas* c1 = new TCanvas("glc1", "TH2Poly from a kdTree",0,0,600,800);
   c1->Divide(1,3);
   c1->cd(1);
   h1->Draw("lego");

   c1->cd(2);
   h2pol->Draw("COLZ L");
   c1->Update();


   /* Draw an equivalent plot showing the data points */
   /*-------------------------------------------------*/

   std::vector<Double_t> z = std::vector<Double_t>(DATASZ, 0.);
   for (UInt_t i = 0; i < DATASZ; ++i)
      z[i] = (Double_t) h2pol->GetBinContent(h2pol->FindBin(smp[i], smp[DATASZ + i]));

   TGraph2D *g = new TGraph2D(DATASZ, smp, &smp[DATASZ], &z[0]);
   gStyle->SetPalette(1);
   g->SetMarkerStyle(20);

   c1->cd(3);
   g->Draw("pcol");
   c1->Update();

   // ---------------------------------------------------------
   // make a new TH2Poly where bins are ordered by the density
   // ---------------------------------------------------------

   TH2Poly* h2polrebin = new TH2Poly("h2PolyBinTest", "KDTree binning", kdBins->GetDataMin(0), kdBins->GetDataMax(0), kdBins->GetDataMin(1), kdBins->GetDataMax(1));
   h2polrebin->SetFloat();

   /*---------------------------------*/
   /* Sort the bins by their density  */
   /*---------------------------------*/

   kdBins->SortBinsByDensity();

   for (UInt_t i = 0; i < kdBins->GetNBins(); ++i) {
      const Double_t* binMinEdges = kdBins->GetBinMinEdges(i);
      const Double_t* binMaxEdges = kdBins->GetBinMaxEdges(i);
      h2polrebin->AddBin(binMinEdges[0], binMinEdges[1], binMaxEdges[0], binMaxEdges[1]);
   }

   for (UInt_t i = 1; i <= kdBins->GetNBins(); ++i){
      h2polrebin->SetBinContent(i, kdBins->GetBinDensity(i - 1));}

   std::cout << "Bin with minimum density: " << kdBins->GetBinMinDensity() << std::endl;
//.........这里部分代码省略.........

void SimRC_eventisuareatotale(int seed)
{
	double x[N_events], y[N_events], z[N_events], zfin[N_events], x1final[N_events], y1final[N_events], xfin[N_events], yfin[N_events], phi[N_events], theta[N_events], x1[N_events], y1[N_events], z1[N_events], phi1[N_events], theta1[N_events]; // all the necessary array
	double uniform, u;
	TRandom3 *rand = new TRandom();
	rand->SetSeed(seed); // to set the seed
	int counterge = 0;		//to count good events
	int k = 0;		//dynamic counter 
	for (int i = 0; i < N_events; i++) {
		//rand->SetSeed(0);		//forget about this, it's wrong; I kept it only to remember not to use it
		x[i] = rand->Uniform(-square,L+square);		//generates N_events values along X
		y[i] = rand->Uniform(-square, D+square);
		phi[i] = rand->Uniform(0, 2*PI);
		uniform = 0;
		u = 2;
		while (M*u > M*cos(uniform)*cos(uniform)*sin(uniform)) {		//rejection condition
			u = rand->Uniform(0, 1);
			uniform = rand->Uniform(0, PI / 2);
		};
		theta[i] = uniform;
		if (x[i] >= 0 && x[i] <= L && y[i] >= 0 && y[i] <= D) {
			z[i] = 0;
			x1[counterge] = x[i];
			y1[counterge] = y[i];
			z1[counterge] = z[i];
			theta1[counterge] = theta[i];
			phi1[counterge] = phi[i];
			x1final[counterge] = x[i];
			y1final[counterge] = y[i];
			counterge++; 
			k++;
			//cout << k << endl;
		}
		else if (x[i] < 0 && phi[i] >PI / 2 && phi[i] < PI && (y[i] - x[i] * tan(phi[i])) > 0 && (y[i] - x[i] * tan(phi[i])) < D && theta[i] > atan(x[i] / (cos(phi[i]) * sp))) {
			z[i] = -x[i] / (cos(phi[i]) * tan(theta[i]));
			x1[counterge] = x[i];
			y1[counterge] = y[i];
			z1[counterge] = z[i];
			theta1[counterge] = theta[i];
			phi1[counterge] = phi[i];
			x1final[counterge] = 0;
			y1final[counterge] = (y[i] - x[i] * tan(phi[i]));
			counterge++;
			k++;
			//cout << k << endl;			
		}
		else if (x[i] < 0 && phi[i] >PI / 2 && phi[i] < PI && (y[i] + x[i] * tan(PI-phi[i])) > D && x[i]-(D-y[i])/ tan(PI-phi[i]) > 0 && x[i] - (D - y[i]) / tan(PI-phi[i]) <L && theta[i] > atan((y[i]-D) / (sin(PI-phi[i]) * sp))) {
			z[i] = -(y[i] - D) / (sin(PI - phi[i]) * tan(theta[i]));
			x1[counterge] = x[i];
			y1[counterge] = y[i];
			z1[counterge] = z[i];
			theta1[counterge] = theta[i];
			phi1[counterge] = phi[i];
			y1final[counterge] = D;
			x1final[counterge] = x[i] - (D - y[i]) / tan(PI - phi[i]);
			counterge++;
			k++;
			//cout << k << endl;
		}
		else if (x[i] < 0 && phi[i] >PI && phi[i] < 3*PI/2 && (y[i] - x[i] * tan(phi[i])) > 0 && (y[i] - x[i] * tan(phi[i])) < D && theta[i] > atan(-x[i] / (cos(PI-phi[i]) * sp))) {
			z[i] = -x[i] / (cos(phi[i]) * tan(theta[i]));
			x1[counterge] = x[i];
			y1[counterge] = y[i];
			z1[counterge] = z[i];
			theta1[counterge] = theta[i];
			phi1[counterge] = phi[i];
			x1final[counterge] = 0;
			y1final[counterge] = (y[i] - x[i] * tan(phi[i]));
			counterge++;
			k++;
			//cout << k << endl;
		}
		else if (x[i] < 0 && phi[i] >PI && phi[i] < 3 * PI / 2 && x[i] - y[i] / tan(phi[i]) > 0 && x[i] - y[i] / tan(phi[i]) < L&& (y[i] - x[i] * tan(phi[i])) < 0 && theta[i] > atan(-y[i] / (sin(phi[i]-PI) * sp))) {
			z[i] = y[i] / (sin(phi[i] - PI) * tan(theta[i]));
			x1[counterge] = x[i];
			y1[counterge] = y[i];
			z1[counterge] = z[i];
			theta1[counterge] = theta[i];
			phi1[counterge] = phi[i];
			y1final[counterge] = 0;
			x1final[counterge] = x[i] - y[i] / tan(phi[i]);
			counterge++;
			k++;
			//cout << k << endl;
		}
		else if (x[i] > L && phi[i] > 0 && phi[i] < PI / 2 && (y[i] - (x[i] - L) * tan(phi[i])) > 0 && (y[i] - (x[i] - L) * tan(phi[i])) < D && theta[i] > atan((x[i] - L) / (cos(phi[i]) * sp))) {
			z[i] = -(x[i] - L) / (cos(phi[i]) * tan(theta[i]));
			x1[counterge] = x[i];
			y1[counterge] = y[i];
			z1[counterge] = z[i];
			theta1[counterge] = theta[i];
			phi1[counterge] = phi[i];
			x1final[counterge] = L;
			y1final[counterge] = y[i] - (x[i] - L) * tan(phi[i]);
			counterge++;
			k++;
			//cout << k << endl;
		}
		else if (x[i] > L && phi[i] > 0 && phi[i] < PI / 2 && (y[i] - (x[i] - L) * tan(phi[i])) > D && x[i] - (y[i] - D) / tan(phi[i]) > 0 && x[i] - (y[i] - D) / tan(phi[i]) < L && theta[i] > atan((y[i]-D) / (sin(phi[i]) * sp))) {
			z[i] = -(y[i] - D) / (sin(phi[i]) * tan(theta[i]));
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	
	
	//______________________________
	//   ======================================================================
	//============================================================================
	//START ANALYSIS==============================================================
	//============================================================================
	//   ======================================================================
	//*
	cout << "\n";
	cout << "<----STARTING ANALYSIS---->" << endl;
  cout << "\n";
	cout << pairTree->GetEntries() << " pairs to analyze" << endl;
	
	
	
	int blueFillNo;
	int yellowFillNo;
	
	int phiSRbin;
	
	
	vector<double> fillAsyms;
	vector<double> fillAsymsE;
	vector<double> fillVec;
	
	TLorentzVector sum;
	TLorentzVector sumY;
	TLorentzVector sumB;
	
	//random number for randomizing spin.
	//set seed to zero to gaurenty unique numbers each time.
	TRandom3 r;
	r.SetSeed(0);
	
	int rand5 = 0;
	int rand6 = 0;
	int rand9 = 0;
	int rand10 = 0;
	
	
	int totalPairsFinal = 0;
	int blueD = 0;
	int blueU = 0;
	int yellowD = 0;
	int yellowU = 0;
	
	int pionStarNumber = 0;
	int runsProcessed = 0;
	
	double currentFillNo = 0;
	double currentRunNo = 0;

	
	for (int iPair = pionStarNumber; iPair < pairTree->GetEntries(); iPair++)
	{
		
		if (iPair%10000 == 0) {cout << "processing pair number " << iPair << endl;}
		//if (iPair == pionStarNumber+300000){break;}
		
		pairTree->GetEntry(iPair);
		
		//if (runsProcessed > 4){break;} 
		
		
		if (pair1->runInfo().beamFillNumber(1) != currentFillNo && currentFillNo != 0 && currentFillNo != 16427)

//.........这里部分代码省略.........
    }
	
    // Setup electron transport
    AvalancheMicroscopic* aval = new AvalancheMicroscopic();
    aval->SetSensor(sensor);
    //aval->EnableAvalancheSizeLimit(1000);
    
    sensor->AddElectrode(fm, "readout");
    sensor->AddElectrode(fm, "ions");
    const double tMin = 0.;
    const double tMax = 75.;
    const double tStep = 0.2;
    const int nTimeBins = int((tMax - tMin)/tStep);
    sensor->SetTimeWindow(0., tStep, nTimeBins);
    aval->EnableSignalCalculation();
    
    ViewSignal* signalView = new ViewSignal();
    signalView->SetSensor(sensor);
    TH1D* h; // tmp storage of timing histogram

    // Setup ion transport
    AvalancheMC* iondrift = new AvalancheMC();
    iondrift->SetSensor(sensor);
    iondrift->EnableSignalCalculation();
    iondrift->SetDistanceSteps(2e-4);

    // Calculate avalanche
    int ne, ni, np, status, nc;
    double ec, extra;
    double x0, y0, z0, t0, e0, x1, y1, z1, t1, e1, x2, y2, z2, t2, e2, x3, y3, z3, t3, e3;
    double vx0, vy0, vz0, vx1, vy1, vz1;
	
    TRandom3 r;
    r.SetSeed(0);

    TString savefile = workingdir + "output/" + output + ".root";
    TFile f(savefile, "recreate");
    TDirectory *dir = f.mkdir("signals");
    dir->cd();
    
    // Prepare tree for charged particle and clusters
    particle p;
    TTree *pTree = new TTree("pTree", "Charged particle");
    pTree->Branch("x0", &p.x0);
    pTree->Branch("y0", &p.y0);
    pTree->Branch("z0", &p.z0);
    pTree->Branch("vx0", &p.vx0);
    pTree->Branch("vy0", &p.vy0);
    pTree->Branch("vz0", &p.vz0);
    pTree->Branch("t0", &p.t0);
    pTree->Branch("e0", &p.e0);
    pTree->Branch("type", "TString", &p.type);
    pTree->Branch("noClusters", &p.noClusters);
    pTree->Branch("stoppingpower", &p.stoppingPower);
    pTree->Branch("lambda", &p.lambda);
    pTree->Branch("clusters", "std::vector<cluster>", &p.clusters);
  
    // Prepare tree for electrons
    avalancheE aE;
    TTree *ETree = new TTree("eTree", "Avalanche electrons");
    ETree->Branch("x0", &aE.x0);
    ETree->Branch("y0", &aE.y0);
    ETree->Branch("z0", &aE.z0);
    ETree->Branch("vx0", &aE.vx0);
    ETree->Branch("vy0", &aE.vy0);
    ETree->Branch("vz0", &aE.vz0);

void estimate_uncertainty(int N=1000)
{
  double a,b,c,d;
  TRandom3 rand;
  double mean = 0.0;
  double sigma = 0.25;
  double sigma1 = 0.22;

  TH1F *hA = new TH1F("hA","numerator 1;A;counts", 200,0,2);
  TH1F *hB = new TH1F("hB","denominator 1;B;counts", 200,0,2);
  TH1F *hC = new TH1F("hC","numerator 2;C;counts", 200,0,2);
  TH1F *hD = new TH1F("hD","denominator 2;D;counts", 200,0,2);
  TH1F *hR1 = new TH1F("hR1","ratio 1;single ratio;counts",200,0,2);
  TH1F *hR2 = new TH1F("hR2","ratio 2;single raito;counts",200,0,2);
  TH1F *hDR = new TH1F("hDR","double ratio;double ratio;counts",200,0,2);

  hA->Sumw2();
  hB->Sumw2();
  hC->Sumw2();
  hD->Sumw2();
  hR1->Sumw2();
  hR2->Sumw2();
  hDR->Sumw2();

  hA->SetMarkerColor(kBlack);
  hB->SetMarkerColor(kRed);
  hC->SetMarkerColor(kBlue);
  hD->SetMarkerColor(kGreen+2);

  hB->SetMarkerStyle(21);
  hC->SetMarkerStyle(24);
  hD->SetMarkerStyle(25);

  hR2->SetMarkerColor(kRed);
  hR2->SetMarkerStyle(24);

  hA->GetXaxis()->CenterTitle(1);
  hR1->GetXaxis()->CenterTitle(1);
  hDR->GetXaxis()->CenterTitle(1);

  rand.SetSeed(0);

  for (int i=0; i<N; ++i) {
    // a = rand.Gaus(mean, sigma1);
    // b = a;
    a = rand.Gaus(mean, sigma1);
    a *= rand.Gaus(mean, sigma);
    b = rand.Gaus(mean, sigma1);
    b *= rand.Gaus(mean, sigma);

    // c = rand.Gaus(mean, sigma1);
    // d = c;
    c = rand.Gaus(mean, sigma1);
    c *= rand.Gaus(mean, sigma);
    d = rand.Gaus(mean, sigma1);
    d *= rand.Gaus(mean, sigma);

    a += 1;
    b += 1;
    c += 1;
    d += 1;

    hA->Fill(a);
    hB->Fill(b);
    hC->Fill(c);
    hD->Fill(d);

    hR1->Fill( (a/b) );
    hR2->Fill( (c/d) );

    hDR->Fill( (a/b) / (c/d) );
  }

  TCanvas *c1 = new TCanvas("c1","c1",1200,600);
  c1->Divide(3,1);
  c1->cd(1);
  hA->Draw();
  hB->Draw("same");
  hC->Draw("same");
  hD->Draw("same");
  c1->cd(2);
  hR1->Draw();
  hR2->Draw("same");
  c1->cd(3);
  hDR->Draw();

  cout << "Mean = " << hDR->GetMean() << " RMS = " << hDR->GetRMS() << endl;

  return;
}

//.........这里部分代码省略.........
		polOfBinSumDown[i] = 0;
		
		avgPerrorOfBinDown[i] = 0;
		pErrorOfBinDown[i] = 0;
		
	}
	
	
	//   ======================================================================
	//============================================================================
	//START ANALYSIS==============================================================
	//============================================================================
	//   ======================================================================
	
	cout << pairTree->GetEntries() << endl;
	
	cout << "\n";
	cout << "<----STARTING ANALYSIS---->" << endl;
  cout << "\n";
	
	
	double blueFillNo;
	double yellowFillNo;
	
	int bin;
	
	TLorentzVector sum;
	TLorentzVector sumY;
	TLorentzVector sumB;
	
	//random number for randomizing spin.
	//set seed to zero to gaurenty unique numbers each time.
	TRandom3 r;
	r.SetSeed(0);
	
	int rand5 = 0;
	int rand6 = 0;
	int rand9 = 0;
	int rand10 = 0;
	
	
	int totalPairsFinal = 0;
	int blueD = 0;
	int blueU = 0;
	int yellowD = 0;
	int yellowU = 0;
	
	int pionStarNumber = 0;
	
	
	for (int iPair = pionStarNumber; iPair < pairTree->GetEntries(); iPair++)
	{
		if (iPair%10000 == 0) {cout << "processing pair number " << iPair << endl;}
		
		//if (iPair == pionStarNumber+2000000){break;}
		
		pairTree->GetEntry(iPair);
				
		
		if (pair1->withinRadius(0.05, 0.3))
		{
			
			
			bool triggerFired = false;
			bool fromKaon = false;
			

//.........这里部分代码省略.........
	
	cout << "Pt  between " << lowLimitPt << " and " << hiLimitPt << endl;
	cout << "M   between " << lowLimitMass << " and " << hiLimitMass << endl;
	cout << "Eta between " << lowLimitEta << " and " << hiLimitEta << endl;
	//______________________________
	
	//*
	
	//   ======================================================================
	//============================================================================
	//START ANALYSIS==============================================================
	//============================================================================
	//   ======================================================================
		//*
	cout << "\n";
	cout << "<----STARTING ANALYSIS---->" << endl;
  cout << "\n";
	cout << pairTree->GetEntries() << " pairs to analyze" << endl;
	
	
	
	int blueFillNo;
	int yellowFillNo;
	
	int phiSRbin;
	
	TLorentzVector sum;
	TLorentzVector sumY;
	TLorentzVector sumB;
	
	//random number for randomizing spin.
	//set seed to zero to gaurenty unique numbers each time.
	TRandom3 r;
	r.SetSeed(0);
	
	int rand5 = 0;
	int rand6 = 0;
	int rand9 = 0;
	int rand10 = 0;
	
	
	int totalPairsFinal = 0;
	int blueD = 0;
	int blueU = 0;
	int yellowD = 0;
	int yellowU = 0;
	
	int pionStarNumber = 0;
	int test = 0;
	
	for (int iPair = pionStarNumber; iPair < pairTree->GetEntries(); iPair++)
	{
		
		if (iPair%10000 == 0) {cout << "processing pair number " << iPair << endl;}
		//if (iPair == pionStarNumber+100000){break;}

		pairTree->GetEntry(iPair);

		
		if (pair1->withinRadius(0.05, 0.3))
		{

			bool triggerFired = false;
			bool fromKaon = false;
			bool passInitialCut_B = false;
			bool passInitialCut_Y = false;