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

void macro6(){

    TH1F* sig_h=new TH1F("sig_h","Signal Histo",50,0,10);
    TH1F* gaus_h1=new TH1F("gaus_h1","Gauss Histo 1",30,0,10);
    TH1F* gaus_h2=new TH1F("gaus_h2","Gauss Histo 2",30,0,10);
    TH1F* bkg_h=new TH1F("exp_h","Exponential Histo",50,0,10);

    // simulate the measurements
    TRandom3 rndgen;
    for (int imeas=0;imeas<4000;imeas++){
        bkg_h->Fill(rndgen.Exp(4));
        if (imeas%4==0) gaus_h1->Fill(rndgen.Gaus(5,2));
        if (imeas%4==0) gaus_h2->Fill(rndgen.Gaus(5,2));
        if (imeas%10==0)sig_h->Fill(rndgen.Gaus(5,.5));}

    // Format Histograms
    TH1F* histos[4]={sig_h,bkg_h,gaus_h1,gaus_h2};
    for (int i=0;i<4;++i){
        histos[i]->Sumw2(); // *Very* Important
        format_h(histos[i],i+1);
        }

    // Sum
    TH1F* sum_h= new TH1F(*bkg_h);
    sum_h->Add(sig_h,1.);
    sum_h->SetTitle("Exponential + Gaussian");
    format_h(sum_h,kBlue);

    TCanvas* c_sum= new TCanvas();
    sum_h->Draw("hist");
    bkg_h->Draw("SameHist");
    sig_h->Draw("SameHist");

    // Divide
    TH1F* dividend=new TH1F(*gaus_h1);
    dividend->Divide(gaus_h2);

    // Graphical Maquillage
    dividend->SetTitle(";X axis;Gaus Histo 1 / Gaus Histo 2");
    format_h(dividend,kOrange);
    gaus_h1->SetTitle(";;Gaus Histo 1 and Gaus Histo 2");
    gStyle->SetOptStat(0);

    TCanvas* c_divide= new TCanvas();
    c_divide->Divide(1,2,0,0);
    c_divide->cd(1);
    c_divide->GetPad(1)->SetRightMargin(.01);
    gaus_h1->DrawNormalized("Hist");
    gaus_h2->DrawNormalized("HistSame");

    c_divide->cd(2);
    dividend->GetYaxis()->SetRangeUser(0,2.49);
    c_divide->GetPad(2)->SetGridy();
    c_divide->GetPad(2)->SetRightMargin(.01);
    dividend->Draw();
}

//.........这里部分代码省略.........
      Int_t index=0;
      MatchObject m;
      vector<MatchObject> matches;
      while (index>=0 && index<pars.max_distance_unit_propagation && index<Max(firstmatchx,firstmatchy)){
	if (propray(lr,kNS,index,matchx,m)){
	  if (!matched) firstmatchx=matchx;
	  matched=true;
	  matches.push_back(m);
	}
	if (propray(lr,kEW,index,matchy,m)){
	  if (!matched) firstmatchy=matchy;
	  matched=true;
	  matches.push_back(m);
	}
	index++;
      }

      vector<GoodMatchObject> goodmatches=convert_matches(matches);
	
      Double_t pathxy = pars.max_distance_unit_propagation*10*pars.xtalsize;
      GoodMatchObject finalmatch;
      for (size_t i=0; i<goodmatches.size(); i++){
	Double_t path = sqrt(pow(goodmatches[i].positionx-lr.origin.x(),2)+pow(goodmatches[i].positiony-lr.origin.y(),2));
	if (path<pathxy) {pathxy=path; finalmatch=goodmatches[i];}
      }
      if (pathxy>1e4){
	if (debug) cout << "LOOPER" << endl;
	continue;
      }

      Int_t channel = findchannel(finalmatch.x,finalmatch.y)-1;
     
      if (debug) cout << finalmatch.x << " " << finalmatch.y << " " << pathxy << " " << channel+1 << endl;

      Double_t path3d = pathxy/Sin(lr.dirvect.Theta());
     
      TVector3 rotated_origin_xy = lr.origin; rotated_origin_xy.SetZ(0);
      TVector3 rotated_impact = TVector3(finalmatch.positionx,finalmatch.positiony,0);
      TRotation r;
      r.RotateZ(-Pi()/4);
      rotated_origin_xy = r*rotated_origin_xy;
      rotated_impact = r*rotated_impact;

      Int_t nx = 0;
      Int_t ny = 0;
      Int_t nz = 0;
      {
	Int_t sx = finalmatch.x+finalmatch.y;
	if (sx==0) nx=0;
	else if (sx>0) nx = sx/2-1;
	else nx = TMath::Abs(sx)/2;
	Int_t dy = finalmatch.y-finalmatch.x;
	if (dy==0) ny=0;
	else if (dy>0) ny = dy/2-1;
	else ny = TMath::Abs(dy)/2;
	Float_t finalz = path3d*Cos(lr.dirvect.Theta())+deposit.position.z()-pars.xtalheight/2;
	nz = Int_t((fabs(finalz)+pars.xtalheight/2)/pars.xtalheight);
      }
      Int_t all_crossings = nx+ny+nz;

      Int_t my_paper_refl = 0;

      TVector2 difference = TVector2(fabs(rotated_origin_xy.x()-rotated_impact.x()),fabs(rotated_origin_xy.y()-rotated_impact.y()));
      if (difference.Phi()<limit_angle) my_paper_refl+=nx;
      if (Pi()/2-difference.Phi()<limit_angle) my_paper_refl+=ny;
      if (lr.dirvect.Theta()<limit_angle) my_paper_refl+=nz;

      Double_t att_absorption_point = randgen.Exp(pars.light_att_length);
      if (path3d>att_absorption_point) continue;

      Double_t prob_loss_in_reflections = pow(pars.eff_reflection_paper,my_paper_refl)*pow(pars.eff_reflection_total,all_crossings-my_paper_refl);
      if (randgen.Uniform()>prob_loss_in_reflections*pars.eff_lightcoll) continue;
	
      Double_t scintillation_delay = randgen.Exp(pars.scintillation_typ_timescale);

      Double_t arrival_time = deposit.time+path3d/pars.speed_of_light+scintillation_delay;
      if (arrival_time>pars.limit_arrival_time) continue;

      myphotons.at(channel)++;
      output->chamfer_pulseshape[channel]->Fill(arrival_time);
      output->reflections_paper->Fill(my_paper_refl);
      output->reflections_total->Fill(all_crossings-my_paper_refl);
      output->reflections_all->Fill(all_crossings);
      output->optical_path->Fill(path3d);
	
    }

    for (int i=0; i<4; i++) output->chamfer_photons[i]->Fill(myphotons[i]);

  }

  for (int i=0; i<4; i++) Normalize(output->chamfer_pulseshape[i]);
  Normalize(output->reflections_paper);
  Normalize(output->reflections_total);
  Normalize(output->reflections_all);
  Normalize(output->optical_path);

  return output;

};