C++ TF1::GetRandom方法代码示例

void test(){
TH1D* h = new TH1D("","",1000,-5,5);
for(int i=0;i<1e5;i++){

   if(i%10000==0) cout<<i<<endl;
  // TF1 * f = new TF1("f","x",0,1);
  TF1 *f = new TF1("f","exp(-(x-2.1)^2/6.3)+exp(-(x+2.1)^2/6.3)",-2.4,2.4);
//TRandom3 *r1 = new TRandom3(0);
//double x = r1->Rndm();
double x = f->GetRandom();
//TRandom3 *r2 = new TRandom3(0);
//double y = r1->Rndm();
double y = f->GetRandom();
h->Fill(x-y);
}
h->Draw();
}

//g++ `root-config --cflags --libs` altgraph.C -o altgraph.out
void altgraph()
{
  TH1F *h1 = new TH1F("h1","gaussian",100,-5,10);
  Float_t sigma,mpv;
  sigma = 3.; mpv = 2.;
  TF1 *fgaus = new TF1("fgaus","exp(-((x-[0])/[1])^2)*(x>0)",-5,10);
  fgaus->SetParameter(0,mpv);
  fgaus->SetParameter(1,sigma);
  h1->FillRandom("fgaus",1000000);
  TCanvas *c1 = new TCanvas();
  h1->Draw();
  for(int p=0; p<20; p++)
  cout<<fgaus->GetRandom()<<endl;
  c1->SaveAs("alt.pdf");
}

void GetRandomTest(){
	double k0=1.39;
	double k1 = 0.425;
	double theta0 = 3.41;
	double theta1 = 1.30;
	int iNpart=2;
	double k_=k0+k1*(iNpart-2);
        double theta_=theta0+theta1*TMath::Log(iNpart-1);
	TF1 *f = new TF1("f","TMath::GammaDist(x,[0],0,[1])",0,200);
	f->SetParameters(k1,theta_);
	cout<<"Value at 0 = "<<f->Eval(0)<<endl;
	cout<<"Value at 1e-11 = "<<f->Eval(1e-11)<<endl;
	cout<<"Integral 1= "<<f->Integral(f->GetXmin(),f->GetXmax())<<endl;
	f->SetRange(1e-12,200);
	cout<<"Integral 2= "<<f->Integral(f->GetXmin(),f->GetXmax())<<endl;
	cout<<"fXmin = "<<f->GetXmin()<<"\tfXmax = "<<f->GetXmax()<<"\tfNpx = "<<f->GetNpx()<<endl;
	f->SetNpx(1e5);
	TCanvas *c1 = new TCanvas();
	c1->SetLogy();
	cout<<"f mean = "<<f->Mean(0,200)<<endl;
	cout<<"math mean = "<<f->GetParameter(0)*f->GetParameter(1)<<endl;
	TH1D* h = new TH1D("h","h",1000,0,200);
	for(int i=0;i<1e6;i++){
		double para = f->GetRandom();
		h->Fill(para);
	}
	h->Scale(1.0/h->Integral()*1000/200);
	h->GetYaxis()->SetRangeUser(1e-10,1);
	h->SetMarkerStyle(24);
	h->SetMarkerColor(4);
	h->SetMarkerSize(1.1);
	TLegend *leg = new TLegend(0.6,0.7,0.8,0.9);
        leg->SetFillColor(0);
        leg->SetFillStyle(0);
        leg->SetBorderSize(0);
        leg->SetTextFont(42);
        leg->SetTextSize(0.03);
        leg->AddEntry(f,"function","lp");
        leg->AddEntry(h,"filled histogram","lp");
	h->Draw("P");
	f->Draw("same");
        leg->Draw("same");
	cout<<"h mean = "<<h->GetMean(1)<<endl;
	c1->Print("TestGetRandom.png");
	}

void Chi(int num, double* arr)
{
  //double chi[3] = arr;
  //double* pointer = χ
  //if(num > -1)
    {
      TF1* f = new TF1("f", "5*(TMath::Power(x,-4))",1,100);
      TH1D* h1 = new TH1D("h1","h1", 100, -0.5, 99.5);

      for(int i = 0; i < 1000000; i++)
	{
	  h1->Fill(f->GetRandom(1,100));
	  //cout << f.GetRandom(1,100) << endl;
	}
      
      //cout << __FILE__ << __LINE__<< endl;

      TF1* pwrLw = new TF1("pwrLw", "[0]*TMath::Power(x,[1])",1,1000);//changed to 0 parameter
      pwrLw->SetParameter(0,5);
      pwrLw->SetParameter(1,-4);
      //cout << pwrLw->GetParameter(0) << "\t" << pwrLw->GetParameter(1) << endl;

      //cout << __FILE__ << __LINE__<< endl;
      
      h1->Fit("pwrLw","0","",2,25);
      TF1* func = h1->GetFunction("pwrLw");
      double chi = func->GetChisquare();
      double dof = func->GetNDF();
      arr[0] = chi/dof;
      double par0 = func->GetParameter(0);
      double par1 = func->GetParameter(1);
      //cout << par0 << "\t" << par1 << endl;
      arr[1] = CalcChiSqr(h1,par0,par1,2,25);
      //chi[2] = CalcChiSqr(h1,5,-4,2,25);
      //cout << "Outputing values from Chi()" << endl;
      //cout << arr[0] << "\t" << arr[1] << endl;//"\t" << chi[2] << endl;
      //DoCleanUp(h1, f, pwrLw);
    }
    //arr = chi;
    //return chi;
}

//.........这里部分代码省略.........
  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){
      printf("ERROR: Number of hits in event is negative!!!\n");
      break;
    }

    //distribute the hits in the event over the roc accorsing to gaus distribution
    /*
    for(int i = 0; i < hitsInEvent; ++i){

      //random row value
      randomrow->SetSeed(0);
      row = randomrow->Gaus(40,10);
      
      //random column value
      randomcol->SetSeed(0);
      col = randomcol->Gaus(25,8);

      //printf("row: %d | col: %d | adc: %d\n",row,col,adc);

      bpixTree[2]->Fill();

    }
    */


    for(int i = 0; i < hitsInEvent; ++i){

      //random row value
      row = fy->GetRandom();
      
      //random column value
      col = fx->GetRandom();

      //random adc value
      adc = fadc->GetRandom();

      //printf("row: %d | col: %d | adc: %d\n",row,col,adc);


      bpixTree[1]->Fill();

    }


    ++ eventNr;
  }


  printf("Total number of Hits: %d\n",hitsTotal);
  printf("Writing output file: %s \n", outputFileName);

  outputFile->cd();
  outputFile->Write();
  outputFile->Close();

  printf("DONE!\n");

}

void proSTEGvnwithnfv3()
{
  
  int mult = atoi(getenv("MULT"));

  int tot_num=50000;  //50k events
  double MeanMult=(double)mult;
  double RMSMult=10;
  
  int ifile = atoi(getenv("IFILE")); 

  //simple toy event generator
  //TFile f(Form("/lio/lfs/cms/store/user/qixu/flow/NewSTEG/pPbDataV205m%d/vndata_50k_%d.root",mult,ifile), "RECREATE","ROOT file with histograms & tree");
  TFile f(Form("/tmp/xuq7/pPbDataV205m%d/vndata_50k_%d.root",mult,ifile), "RECREATE","ROOT file with histograms & tree");
//  TFile f(Form("vndata_50k_%d.root",mult,ifile), "RECREATE","ROOT file with histograms & tree");
  TTree *tree = new TTree("tree","Event tree with a few branches");
//  tree->Branch("npg", &b_npg, "npg/I");   // # of particles;
//  tree->Branch("phirg", &b_phirg, "phirg/F");  // RP angle;
  tree->Branch("n", &b_n, "n/I");          // same as npg;
  tree->Branch("ptg", &b_ptg, "ptg[n]/F");  // ;
  tree->Branch("etag", &b_etag, "etag[n]/F");
  tree->Branch("phig", &b_phig, "phig[n]/F");
  
  TF1 *EtaDistr = new TF1("EtaDistr","exp(-(x-2.1)^2/6.3)+exp(-(x+2.1)^2/6.3)",-4,4);
  //TF1 *PhiDistr = new TF1("PhiDistr","1+2*[0]*cos(x)+2*[1]*cos(2*x)+2*[2]*cos(3*x)+2*[3]*cos(4*x)+2*[4]*cos(5*x)+2*[5]*cos(6*x)",0,2.*TMath::Pi());
  TF1 *PhiDistr = new TF1("PhiDistr","1+2*[0]*cos(2*x)+2*[1]*cos(3*x)+2*[2]*cos(4*x)+2*[3]*cos(5*x)+2*[4]*cos(6*x)",0,2.*TMath::Pi());
  //TF1 *PtDistr  = new TF1("PtDistr","exp (-(x/.40))+0.0015*exp (-(x/1.5))", 0.2,10);	//V~0.12
  //TF1 *PtDistr  = new TF1("PtDistr","exp (-(x/0.90))+0.15*exp (-(x/15))", 0.1,10);	//V~=0.06
  TF1 *PtDistr  = new TF1("PtDistr","0.03*(exp (-(x/0.594540))+0.00499506*exp (-(x/1.89391)))", 0.3,6.0);	//Real Data
  //  TF1 *PtDistr = new TF1("PtDistr","[0]*x*TMath::Power(1+(sqrt(x*x+[1]*[1])-[1]/[2]),-[3])",0.2,10);
	//PtDistr->SetParameters(118.836,-0.335972,0.759243,118.836);	//Real data fit with Tsallis
  //TF1 *V1vsEta = new TF1("V1vsEta","-exp(-(x+1)^2)/20-x/30+exp(-(x-1)^2)/20",-2.4,2.4); 
  //TF1 *V2vsPt   = new TF1("V2vsPt","((x/3)^1.8/(1+(x/3)^1.8))*(.00005+(1/x)^0.8)",0.2,10);
  //TF1 *V2vsPt = new TF1("V2vsPt","((x/[0])^[1]/(1+(x/[2])^[3]))*(.00005+(1/x)^[4])",0.1,10);
 //	V2vsPt->SetParameters(4.81159,1.80783,3.69272,3.11889,0.931485);	//Real data V~0.05
  //	V2vsPt->SetParameters(5,1.8,3,1.8,0.8); //V~0.06
  TF1 *V2vsPt = new TF1("V2vsPt","((x/3.31699)^2.35142/(1+(x/3.49188)^3.54429))*(.00005+(1/x)^1.50600)",0.3,6.0);
  TF1 *V3vsPt = new TF1("V3vsPt","((x/3.2)^2.3/(1+(x/3.4)^2.1))*(.00005+(1/x)^1.4)",0.3,6.0);
  TF1 *V4vsPt = new TF1("V4vsPt","((x/4.8)^2.1/(1+(x/3.4)^2.1))*(.00005+(1/x)^1.4)",0.3,6.0);
  TF1 *V5vsPt = new TF1("V5vsPt","((x/6.0)^3.2/(1+(x/11.4)^2.1))*(.00005+(1/x)^1.4)",0.3,6.0);
  TF1 *V6vsPt = new TF1("V6vsPt","((x/5.6)^2.4/(1+(x/4.7)^2.1))*(.00005+(1/x)^1.4)",0.3,6.0);
 
    UInt_t iniseed = SetSeedOwn();
    gRandom->SetSeed(iniseed);
  TRandom3 *rnd = new TRandom3(0);
  
  double v1, v2, v3, v4, v5, v6, ph, myphi, mypt, myeta, phi0, Psi;
  int nnf,k;
  double neta, nphi, npt;
  int slicept;
  
  for(int i=0; i<tot_num; i++){ 
    
    Psi = rnd->Uniform(0.0,2.*TMath::Pi());
    //Psi=0;
    b_phirg = Psi; 
    b_npg = rnd->Gaus(MeanMult,RMSMult); 
    int b_npgsame = (int)b_npg * 0.10;
    n = 0;
  
    for(int j=0; j<b_npg;j++ ){
      mypt = PtDistr->GetRandom();
      myeta =  EtaDistr->GetRandom();

      //v1=V1vsEta->Eval(myeta);
      v2=V2vsPt->Eval(mypt);
      v3=V3vsPt->Eval(mypt);
      v4=V4vsPt->Eval(mypt);
      v5=V5vsPt->Eval(mypt);
      v6=V6vsPt->Eval(mypt);

      b_etag[n] = myeta;
      b_ptg[n]  = mypt;
      //PhiDistr->SetParameters(v1,v2,v3,v4,v5,v6);
      PhiDistr->SetParameters(v2,v3,v4,v5,v6);
      
      myphi = PhiDistr->GetRandom(); // randon selection dn/dphi

      myphi = myphi+Psi; // angle in lab frame -- needed for correct cumulant v2
      if (myphi>2.*TMath::Pi()) myphi=myphi-2.*TMath::Pi(); // 0 - 2*Pi
      
      b_phig[n] = myphi; // save angle in lab frame
      n++;
    }
    if(i%10==0){
    for(int j=0;j<b_npgsame;j++){
      mypt = PtDistr->GetRandom();
      myeta =  EtaDistr->GetRandom();
      b_ptg[n] = mypt;
      b_etag[n] = myeta;
      myphi = rnd->Gaus(Psi, 0.3);
      if (myphi>2.*TMath::Pi()) myphi=myphi-2.*TMath::Pi(); // 0 - 2*Pi
      b_phig[n] = myphi;
      n++;
    }
    }
    if (i%10000 == 0) cout << i << " " << "events processed" << endl;

    b_n = n;
    tree->Fill();
//.........这里部分代码省略.........

// Evaluates a transfer function attached to an object
// the tf is a TF1* in a TFType->TF1* map in the object
// the tf is a function of the reconstructed Energy (pt) through tf->Eval(Erec)
// The generator energy is set via tf->SetParameter(0, Egen)
// type - the hypothesis for the object to be tested, e.g. qReco (reconstructed light quark)
double MEM::transfer_function2( void* obj, //either MEM::Object or TF1
				const double *x,
				const TFType::TFType& type,
				int& out_of_range,
				const double& cutoff,
				const bool& smear,
				const int& debug){
  
  double w = 1.0;
  
  TF1* tf = nullptr;
  MEM::Object* _obj = nullptr;
  
  //x[0] -> Egen
  switch( type ){
    
    //W(Erec | Egen) = TF1(Erec, par0:Egen)
  case TFType::bReco:
  case TFType::qReco:
    _obj = (MEM::Object*)obj;
    tf = _obj->getTransferFunction(type);
    
    //set gen
    tf->SetParameter(0, x[0]);
    //eval with x - reco
    w *= tf->Eval(_obj->p4().Pt());
    if( smear )
      w = tf->GetRandom();

#ifdef DEBUG_MODE
    if( debug&DebugVerbosity::integration) 
      cout << "\t\ttransfer_function2: Evaluate W(" << x[0] << " | " << _obj->p4().Pt() << ", TFType::qReco) = " << w << endl;
#endif
    break;
    
    //In case jets have a cut E > Emin
    //W_loss(Egen) = \int_0^Emin W(Erec | Egen) dErec
    //In case jets have a cut pt > ptmin
    //W_loss(ptgen) = \int_0^ptmin W(ptrec | ptgen) dptrec
    //FIXME pt to be implemented
  case TFType::bLost:
  case TFType::qLost:
    
    tf = (TF1*)obj;
    
    //eval at x - gen
    w *= tf->Eval(x[0]);
#ifdef DEBUG_MODE
    if( debug&DebugVerbosity::integration) 
      cout << "\t\ttransfer_function2: Evaluate W(" << x[0] 
	   << ", TFType::qLost) = " << w << endl;
#endif
    break;
    
  case TFType::Unknown:
    w *= 1.;
#ifdef DEBUG_MODE
    if( debug&DebugVerbosity::integration) 
      cout << "\t\ttransfer_function2: Evaluate W = 1 " << endl;
#endif
    break;

  default:
    break;
  }

  return w;
}

void toymc(double ptmin = 290., double ptmax = 320., double lumi = 0,
	   double *nevts = 0, double *mean = 0, double *err = 0) {

  gStyle->SetOptStat(1110);
  gStyle->SetOptFit();

  TF1 *cb1 = new TF1("cb1",crystal_ball, 0., 2., 5);
  TF1 *cb1x = new TF1("cb1x",crystal_ball_x, 0., 2., 5);
  TF1 *g = new TF1("g","gaus",0.,2.);
  TF1 *cb2 = new TF1("cb2",crystal_ball, 0., 2., 5);
  TF1 *cb2x = new TF1("cb2x",crystal_ball_x, 0., 2., 5);
  TF1 *cb3 = new TF1("cb3",crystal_ball, 0., 2., 5);
  TF1 *cb3x = new TF1("cb3x",crystal_ball_x, 0., 2., 5);
  TF1 *cb4 = new TF1("cb4",crystal_dipole, 0., 2., 5);
  TF1 *cb4x = new TF1("cb4x",crystal_dipole_x, 0., 2., 5);

  double mu = 0.70; // mu1*mu2
  double mu1 = 0.95;
  double mu2 = mu/0.95;
  double mu3 = 1.00;
  double mu4 = 1.00;
  // Sigmas are relative widths so remember to multiply by mean
  double sig = 0.12; // sqrt(sig1**2 + sig2**2)
  double sig1 = 0.06;
  double sig2 = sqrt(sig*sig-sig1*sig1);
  double sig3 = 0.025;
  double sig4 = 0.030;
  double a1 = 0.7;//1.2;
  double n1 = 15.;
  double a2 = 1.7;
  double n2 = 15.;
  double a3 = 1.0;
  double n3 = 25.;
  double a4 = 0.8;
  double n4 = 35.;

  const bool photbin = false;//true; // Bin in calo photon pT, not parton pT
  double binmin = ptmin;//290;
  double binmax = ptmax;//320;
  double genmin = binmin*0.93;//0.9;
  double genmax = binmax*1.40;//1.20;

  //TF1 *xsec = new TF1("xsec","7.729e+10*pow(x,-3.104)*exp(-0.003991*x)",
  TF1 *xsec = new TF1("xsec", cross_section, 0., 14000., 5);
  xsec->SetParameters(7.729e+10, -3.104, -0.003991, genmin, genmax);
  xsec->SetNpx(1000);

  // Estimate average pT and number of events in bin
  TF1 *xsec_x = new TF1("xsec_x", cross_section_x, 0., 14000., 5);
  xsec_x->SetParameters(7.729e+10, -3.104, -0.003991, genmin, genmax);
  xsec_x->SetNpx(1000);
  double pt = xsec_x->Integral(binmin,binmax) / xsec->Integral(binmin,binmax);
  double nev = xsec->Integral(binmin, binmax) * lumi;
  if (nev==0) {
    nev = 5000000;
    lumi = nev / xsec->Integral(binmin, binmax);
  }
  cout << Form("Processing bin %1.0f-%1.0f with average pT=%1.0f"
	       " and lumi=%1.1f pb-1 (nev=%1.3g)",
	       binmin, binmax, pt, lumi, nev) << endl;

  cb1->SetParameters(1., mu1, sig1*mu1, a1, n1); // Gen/Part
  cb1x->SetParameters(1., mu1, sig1*mu1, a1, n1); // Gen/Part
  //g->SetParameters(1., mu2, sig2*mu2); // Calo/Gen
  cb2->SetParameters(1., mu2, sig2*mu2, a2, n2); // Calo/Part
  cb2x->SetParameters(1., mu2, sig2*mu2, a2, n2); // Calo/Part
  cb3->SetParameters(1., mu3, sig3*mu3, a3, n3); // Calo/Phot
  cb3x->SetParameters(1., mu3, sig3*mu3, a3, n3); // Calo/Phot
  cb4->SetParameters(1., mu4, sig4*mu4, a4, n4); // Part/Part
  cb4x->SetParameters(1., mu4, sig4*mu4, a4, n4); // Part/Part

  double mu1x = cb1x->Integral(0.,2.) / cb1->Integral(0.,2.);
  double mu2x = cb2x->Integral(0.,2.) / cb2->Integral(0.,2.);
  double mu3x = cb3x->Integral(0.,2.) / cb3->Integral(0.,2.);
  double mu4x = cb4x->Integral(0.,2.) / cb4->Integral(0.,2.);

  TH1D *xsec0 = new TH1D("xsec0","PtHat",1000,0.7*genmin,1.3*genmax);
  TH1D *xsec1 = new TH1D("xsec1","PartPtPho",1000,0.7*genmin,1.3*genmax);
  TH1D *xsec2 = new TH1D("xsec2","PhotPt",1000,0.7*genmin,1.3*genmax);
  TH1D *xsec3 = new TH1D("xsec3","PhotPt",1000,0.7*genmin,1.3*genmax);
  TH1D *h1 = new TH1D("h1","GenPart",400,0.,2.);
  TH1D *h2 = new TH1D("h2","CaloGen",400,0.,2.);
  TH1D *h3 = new TH1D("h3","Phot",400,0.,2.);
  TH1D *h4 = new TH1D("h4","Part",400,0.,2.);
  TH1D *h12 = new TH1D("h12","CaloPart",400,0.,2.);
  TH1D *h123 = new TH1D("h123","CaloPhot",400,0.,2.);
  TH1D *h1234 = new TH1D("h1234","RmeasSig",400,0.,2.);

  const double pthat0 = 300.;
  const int nit = 5000000;
  for (int i = 0; i != nit; ++i) {

    double pthat = pthat0;
    //if (photbin) 
    pthat = xsec->GetRandom(genmin, genmax);

    // Is the parton balancing model ok?
    // Should we rather simulate a "second jet", which either
    // the photon or jet can loose? (No gains allowed vs pthat)

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

//.........这里部分代码省略.........
         jetAbove120 = true ;
       }
      if(!jetAbove40)                 Nevents[curr_bin]++ ;
      else if(jetAbove40 && !jetAbove60) Nevt_40_60[curr_bin]++;
      else if(jetAbove60 && !jetAbove75) Nevt_60_75[curr_bin]++;
      else if(jetAbove75 && !jetAbove95) Nevt_75_95[curr_bin]++;
      else if(jetAbove95 && !jetAbove120) Nevt_95_120[curr_bin]++;
     // if(jetAbove120)                  Nevt_120[curr_bin]++;
      else                               Nevt_120[curr_bin]++;
    
      if(TrigName=="Jet20") jetAbove = jetAbove40 && !jetAbove60 ;
      else if(TrigName=="Jet40") jetAbove = jetAbove60 && !jetAbove75 ;
      else   if(TrigName=="Jet60") jetAbove = jetAbove75 && !jetAbove95  ;
      else   if(TrigName=="Jet80") jetAbove = jetAbove95 && !jetAbove120  ;
      else   if(TrigName=="Jet100") jetAbove = jetAbove120  ;
      else jetAbove = !jetAbove40 ;

      if(!jetAbove) continue ;
    
   // for inclusive jet analysis
   for(int j4i = 0; j4i < my_ct->nref ; j4i++) {
    double jetweight = 1; 
    double jet_pt= my_ct->jtpt[j4i];
    double raw_pt= my_ct->rawpt[j4i];
    double jet_eta = my_ct->jteta[j4i];
    if (my_ct->rawpt[j4i]<15) continue;
    int dEtaBin = -1. ;
 //   if( my_hists->IsMC==kTRUE && my_ct->subid[j4i] != 0) continue;
     //for jet kinematcis cuts
     if(TMath::Abs(jet_eta)<=3.){
          my_hists->jetptJES[curr_bin]->Fill(jet_pt, raw_pt/jet_pt, weight);
       if( my_hists->IsMC!=kTRUE ) jetweight*=(fUE->Eval(jet_pt))*C_rel->GetBinContent(C_rel->FindBin(jet_eta));
       else
          jetweight*=((fUE->Eval(jet_pt))*fgaus->GetRandom()+1);
          my_hists->jetptEta[curr_bin]->Fill(jet_pt*jetweight, jet_eta, weight);
   if(coll=="PPb"){       
       if(TMath::Abs(jet_eta+0.465)<=1.) my_hists->jetpt[curr_bin]->Fill(jet_pt*jetweight, weight);
     }
   if(coll=="PbP"){  
       if(TMath::Abs(jet_eta-0.465)<=1.) my_hists->jetpt[curr_bin]->Fill(jet_pt*jetweight, weight);
     } 
       if((jet_pt*jetweight)>100.) my_hists->jetEta[curr_bin]->Fill(jet_eta, weight);
         for(Int_t ieta = 0 ; ieta <netabin; ieta++){
            if(coll=="PPb"){
               if((jet_eta+0.465)>deta[ieta]&&(jet_eta+0.465)<=deta[ieta+1]) dEtaBin = ieta ;
             }
            else if(coll=="PbP"){
              if((jet_eta-0.465)>deta[ieta]&&(jet_eta-0.465)<=deta[ieta+1]) dEtaBin = ieta ;
             }
           else {
             if((jet_eta+0.465)>deta[ieta]&&(jet_eta+0.465)<=deta[ieta+1]) dEtaBin = ieta ;
            } 
         } //assign the eta bin for jets
      if(dEtaBin!=-1){
        my_hists->jetptEtaBin[curr_bin][dEtaBin]->Fill(jet_pt*jetweight, weight);   
        my_hists->NjetsEtaBin[curr_bin][dEtaBin]->Fill(1);
     } 
    }// for jet kinematics cuts
   } //! jet loop

/*        //Leading Jets seach, for tracking efficiency
      double leadingJetPt = -1. ;
      Int_t leadingJetIndex = -1 ;
     for(int j = 0; j < my_ct->nref ; j++) {
         if (fabs(my_ct->jteta[j])>2.5) continue;
         if (my_ct->rawpt[j]<15) continue;

TestProblem
AtlasDiJetMass(const int Nt,
               const int Nr,
               double tbins[],
               double rbins[],
               const double apar,
               const double bpar,
               const int nEvents,
               const double evtWeight)
{
  // From "Fully Bayesian Unfolding" by G. Choudalakis.
  // see arXiv:1201.4612v4
  // Recommended binning:
  // double bins[Nt+1] = {0};
  // for (int j=0; j<=Nt; j++)
  //   bins[j] = 500*TMath::Exp(0.15*j);

  static int id = 0; id++;
  TestProblem t;
  TRandom3 ran;

  // Mass distribution dN/dM
  TF1 *mass = new TF1("mass_dist",
                      "TMath::Power(1.-x/7000,6.0)/TMath::Power(x/7000,4.8)",
                      tbins[0], tbins[Nt]);

  // Generate test problem by MC convolution
  TH1D *hT   = new TH1D("hT",   "Truth mass dist. #hat{T}", Nt, tbins);
  TH1D *hTmc = new TH1D("hTmc", "MC Truth mass dist. #tilde{T}", Nt, tbins);
  TH1D *hD   = new TH1D("hD", "Measured mass dist.", Nr, rbins);
  TH2D *hM   = new TH2D("hM", "Migration matrix", Nr, rbins, Nt, tbins);
  hM->SetTitle(Form("%d x %d migration matrix M_{tr};"
                    "mass (observed);mass (true)", Nr, Nt));
  hT->SetLineWidth(2);
  hD->SetLineWidth(2);

  std::cout << Form("Generating test problem...") << std::flush;
  // Fill histos with MC events.
  // The response matrix gets more statistics than the data.
  double xt, xm, sigma;
  for (int i=0; i<nEvents; i++)
  {
    xt = mass->GetRandom();
    sigma = apar*TMath::Sqrt(xt) + bpar*xt;
    xm = xt + ran.Gaus(0, sigma);
    hM->Fill(xm, xt);
    hTmc->Fill(xt, evtWeight);
    hD->Fill(xm, evtWeight);
  }

  // Simulate Poisson fluctuations in real data (integer content, empty bins)
  for (int r=1; r<=Nr; r++)
    hD->SetBinContent(r, ran.Poisson(hD->GetBinContent(r)));

  // The true truth \hat{T}
  double totmass = mass->Integral(tbins[0],tbins[Nt]);
  for (int j=1; j<=Nt; j++)
  {
    hT->SetBinContent(j, evtWeight*nEvents*mass->Integral(tbins[j-1],
                      tbins[j])/totmass);
    hT->SetBinError(j, 0);
  }
  cout << "Done." << endl;

  // Projection of migration matrix to truth axis. hMt is the
  // numerator for efficiency. Bin contents should be counts
  // here. Elements are normalized to contain probabilities only after
  // projection & division by T-tilde.
  TH1D *hMt  = hM->ProjectionY("hMt",1,Nr);
  TH1D *heff = (TH1D *)hMt->Clone("heff");
  heff->Divide(hTmc);
  heff->Scale(evtWeight);
  hM->Scale(1./nEvents);
  hMt->Scale(1./nEvents);

  t.Response  = hM;
  t.xTruth    = hT;
  t.xTruthEst = hTmc;
  t.xIni      = hMt;
  t.bNoisy    = hD;
  t.eff       = heff;

  return t;
}

void reingold_ROOThomework(){

// Initalizing the Canvas

	TCanvas *c1 = new TCanvas("c1" , "Homework Plots" , 800 , 600 );
	c1 -> Divide(3,3);

// Defining the function with three Gaussians on a quadratic background.

	TF1 *func = new TF1("func", "pol2(0) + gaus(3) + gaus(6) + gaus(9)" , 0 , 100 );
	Double_t param[12] = {9 , 1 , -0.01 , 100 , 20 , 1 , 100 , 50 , 1 , 100 , 80 , 1};
	func->SetParameters(param);
			//  quad, lin, const, A, mu, std , ... 
// Generting and filling the histograms

	TH1F *h0 = new TH1F("h0" , "Original Spectrum" , 100 , 0 , 100 );
	for ( Int_t i = 0 ; i < 2000 ; i++ ){
		Float_t rand = func->GetRandom();
		h0->Fill(rand);
	}
	
	for ( Int_t j = 1 ; j < 4 ; j++ ){
		c1->cd(j);
		h0->Draw();
	}
// Fitting the background using ShowBackground()
	c1->cd(4);
	TH1 *hPeaks = (TH1*) h0->Clone();
	TH1 *hBkgrd = (TH1*) h0->Clone();

	hPeaks->SetName("hPeaks");
	hBkgrd->SetName("hBkgrd");

	hBkgrd = h0->ShowBackground(15);

	h0->Draw();
	hBkgrd->Draw("same");
	
	hPeaks->Add(hBkgrd,-1);
	c1->cd(7);
	hPeaks->Draw();

// Fitting the background without excluding the peaks

	c1->cd(5);
	TH1F *hPeaks2 = (TH1F*) h0->Clone();
	TH1F *hBkg2 = (TH1F*) h0->Clone();

	hPeaks2->SetName("hPeaks2");
	hBkg2->SetName("hBkg2");

	TF1 *quadBack = new TF1("quadBack" , "pol2(0)" , 0 , 100 );
	hBkg2->Fit(quadBack,"R+");
	
	hBkg2->Draw();

	c1->cd(8);

	hPeaks2->Add(quadBack , -1 );
	hPeaks2->Draw();

// Fitting the background excluding the peaks

	c1->cd(6);
	TH1F *hPeaks3 = (TH1F*) h0->Clone();
	TH1F *hBkg3 = (TH1F*) h0->Clone();

	hPeaks3->SetName("hPeaks3");
	hBkg3->SetName("hBkg3");

	TF1 *quadBack2 = new TF1("quadBack2" , fQuad , 0 , 100 ,3);
	hBkg3->Fit(quadBack2,"R+");
	
	hBkg3->Draw();

	c1->cd(9);

	hPeaks3->Add(quadBack2 , -1 );
	hPeaks3->Draw();

}

//.........这里部分代码省略.........
  leg->Draw("same");

  TLatex Tl;
  Tl.SetNDC();
  Tl.SetTextAlign(12);
  Tl.SetTextSize(0.04);
  Tl.SetTextFont(42);
  Tl.DrawLatex(0.18,0.93, "#scale[1.25]{CMS} Performance");
  Tl.DrawLatex(0.65,0.93, Form("%s #sqrt{s_{NN}} = 5.02 TeV",collisionsystem.Data()));

  TLatex* tex;

  tex = new TLatex(0.22,0.78,Form("%.1f < p_{T} < %.1f GeV/c",ptmin,ptmax));
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.04);
  tex->SetLineWidth(2);
  tex->Draw();

  tex = new TLatex(0.22,0.83,"|y| < 1.0");
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.04);
  tex->SetLineWidth(2);
  tex->Draw();

  h->GetFunction(Form("f%d",count))->Delete();
  TH1F* histo_copy_nofitfun = ( TH1F * ) h->Clone("histo_copy_nofitfun");
  histo_copy_nofitfun->Draw("esame");

 TH1D* hTest = new TH1D("hTest","",nbinsmasshisto,minhisto,maxhisto);

 for (int m=0;m<yieldtotal;m++){
 double r = f->GetRandom();
 hTest->Fill(r);
 
 }

 TF1* ffaketotal=(TF1*)f->Clone("ffake");
 TF1* ffakemass=(TF1*)mass->Clone("ffakemass");
 TF1* ffakebackground=(TF1*)background->Clone("ffakebackground");
 TF1* ffakemassSwap=(TF1*)massSwap->Clone("ffakemassSwap");
 
 Double_t yieldtotal_original = ffaketotal->Integral(minhisto,maxhisto)/binwidthmass;
 Double_t yieldmass_original = ffakemass->Integral(minhisto,maxhisto)/binwidthmass;
 Double_t yieldbackground_original = ffakebackground->Integral(minhisto,maxhisto)/binwidthmass;
 Double_t yieldswapped_original = ffakemassSwap->Integral(minhisto,maxhisto)/binwidthmass;

 TH1D* hTestFake = new TH1D("hTestFake","",nbinsmasshisto,minhisto,maxhisto);
  ffakemass->SetParameter(2,ffaketotal->GetParameter(2)*0.8);
  ffakemass->SetParameter(10,ffaketotal->GetParameter(10)*0.8);
 
Double_t yieldmass_modified= ffakemass->Integral(minhisto,maxhisto)/binwidthmass;

cout<<"mass original="<<yieldmass_original<<endl;
cout<<"mass modified="<<yieldmass_modified<<endl;

 for (int m=0;m<yieldmass_original*scalefactor;m++){
   double r = ffakemass->GetRandom();
   hTestFake->Fill(r);
 }

  for (int m=0;m<(int)(yieldbackground_original*scalefactor*bkgreduction);m++){
   double r = ffakebackground->GetRandom();
   hTestFake->Fill(r);
 }

//.........这里部分代码省略.........
//        TCanvas* myCanvas = new TCanvas("myCanvas","myCanvas");
//        htemp->Draw();
//        fsum->Draw("same");

       
//        // loops over toys
       for(int iexp=0; iexp<NEXP; iexp++){

	 TH1D* htoyMC_data = (TH1D*)hdata_data[ieta][ipt]->Clone("htoyMC_data");
	 htoyMC_data->Reset();
	 
	 TH1D* htoyMC_sig  = (TH1D*)hTemplate_S[ieta][ipt]->Clone("htoyMC_sig");
// 	 htoyMC_sig->Reset();
	 
	 TH1D* htoyMC_bkg  = (TH1D*)hTemplate_B[ieta][ipt]->Clone("htoyMC_bkg");
// 	 htoyMC_bkg->Reset();




	 UInt_t nowSeed = (unsigned long)gSystem->Now();
	 gRandom->SetSeed(nowSeed);
	 int nsiggen  = gRandom->Poisson(nsig_input);
	 int nbkggen  = gRandom->Poisson(nbkg_input);
	 int ndata = nsiggen + nbkggen;

	 // reset toy MC data
	 htoyMC_data->Reset();
	 ofstream fout;
	 std::string toyData = Form("/afs/cern.ch/user/s/syu/scratch0/sigSys/toy_%d_%d.dat",ieta,ipt);
	 fout.open(toyData.data());
	 for(int ieve=0; ieve < nsiggen; ieve++)
	   {
	     Double_t xvalue = fsig->GetRandom(fit_lo,fit_hi);
	     fout << xvalue << " " << 
	       0.5*(fBinsPt[ipt]+fBinsPt[ipt+1]) << " " << fBinsEta[ieta] << endl;
	     htoyMC_data->Fill(xvalue);
	   }
	     
	 for(int ieve=0; ieve < nbkggen; ieve++)
	   {
	     Double_t xvalue = fbkg->GetRandom(fit_lo,fit_hi);
	     fout << xvalue << " " << 
	       0.5*(fBinsPt[ipt]+fBinsPt[ipt+1]) << " " << fBinsEta[ieta] << endl;
	     htoyMC_data->Fill(xvalue);
	   }
	     
	 fout.close();

// 	 htoyMC_data->FillRandom("fsig",nsiggen);
// 	 htoyMC_data->FillRandom("fbkg",nbkggen);


// 	 htoyMC_data->Draw();
	 
// 	 cout << "Generated ndata = " << ndata << "\t nsiggen = " << nsiggen << " \t nbkggen = " << 
// 	   nbkggen << endl;


	 Double_t* toyFitResult;
	 Double_t toyMyFitPar[NRETURN]={0};

	 toyFitResult =  Ifit(toyData.data(), 
			      htoyMC_sig, htoyMC_bkg,
			      htoyMC_data, toyMyFitPar,
			      (int)fBinsPt[ipt], dec[ieta],2);

void process_event()
{  
	bool WithProj = true;

	double x_cm = 0;
	double y_cm = 0;

	//Compute centroid
	for(unsigned int j=0; j<particles_targ.size(); j++)
	{
		x_cm += particles_targ[j].x;
		y_cm += particles_targ[j].y;
	}

	double Count = (double) particles_targ.size();

	if (WithProj) {
		for(unsigned int j=0; j<particles_proj.size(); j++)
		{
			x_cm += particles_proj[j].x;
			y_cm += particles_proj[j].y;
		}
		Count += (double) particles_proj.size();
	}

	x_cm = x_cm/Count;
	y_cm = y_cm/Count;

	//Shift origin to centroid
	for(unsigned int j=0; j<particles_targ.size(); j++)
	{
		particles_targ[j].x = particles_targ[j].x - x_cm;
		particles_targ[j].y = particles_targ[j].y - y_cm;
	}

	if (WithProj) {
		//Shift origin to centroid
		for(unsigned int j=0; j<particles_proj.size(); j++)
		{
			particles_proj[j].x = particles_proj[j].x - x_cm;
			particles_proj[j].y = particles_proj[j].y - y_cm;
		}
	}

	//Compute eccentricity
	//JLN - Option to do with blurred Gaussian distributions... (does not change above x_cm,y_cm calculation)
	bool NucleonSmear = true;
	double SmearSigma = 0.4; // units of fm
	TF1 *fradius = new TF1("fradius","x*TMath::Exp(-x*x/(2*[0]*[0]))",0.0,2.0);
	fradius->SetParameter(0,SmearSigma);
	TF1 *fphi = new TF1("fphi","1.0",0.0,2.0*TMath::Pi());
	int SmearSamplings = 100;
	if (NucleonSmear) Count = 0;

	double qx_2 = 0;
	double qy_2 = 0;
	double qx_3 = 0;
	double qy_3 = 0;
	double phi = 0;
	double rsq = 0;
	double r = 0;

	double e2 = 0;
	double e3 = 0;

	//Loop over all particles in the target nucleus and compute eccentricity
	for(unsigned int i=0; i<particles_targ.size(); i++)
	{
		if (!NucleonSmear) {
			r = get_r(particles_targ[i].x,particles_targ[i].y);
			phi = get_phi(particles_targ[i].x,particles_targ[i].y);

			qx_2 += r*r*TMath::Cos(2*phi);
			qy_2 += r*r*TMath::Sin(2*phi);
			qx_3 += r*r*TMath::Cos(3*phi);
			qy_3 += r*r*TMath::Sin(3*phi);
			rsq += r*r;
		}
		if (NucleonSmear) {
			for (int is=0;is<SmearSamplings;is++) {
				double rtemp   = fradius->GetRandom();
				double phitemp = fphi->GetRandom();
				double xtemp   = particles_targ[i].x + rtemp*TMath::Sin(phitemp);
				double ytemp   = particles_targ[i].y + rtemp*TMath::Cos(phitemp);

				r    = TMath::Sqrt(xtemp*xtemp + ytemp*ytemp);
				phi  = TMath::ATan2(ytemp,xtemp);

				qx_2 += r*r*TMath::Cos(2*phi);
				qy_2 += r*r*TMath::Sin(2*phi);
				qx_3 += r*r*TMath::Cos(3*phi);
				qy_3 += r*r*TMath::Sin(3*phi);
				rsq += r*r;

				Count = Count + 1;  // noting this is a double
			}
		}
	}

	if (WithProj) {
//.........这里部分代码省略.........

void Fragment(Int_t nev=0, Int_t debug=0)
{
 gROOT->Reset();
 
 Float_t b;
 Int_t nspectp, nspectn, nspectpfree, nspectnfree; 
 Int_t zz[100], nn[100], nalpha, NFrag, ztot, ntot, ndeu;
 void spectator(Float_t, Int_t*, Int_t*);
 
 TH2F *hsp = new TH2F("hsp","Spectator protons vs b",50,0.,20.,90,0.,90.);
 hsp -> SetXTitle("b (fm)");
 hsp -> SetYTitle("Num. of spectator protons");
 
 TH2F *hsn = new TH2F("hsn","Spectator neutrons vs b",50,0.,20.,130,0.,130.);
 hsn -> SetXTitle("b (fm)");
 hsn -> SetYTitle("Num. of spectator neutrons");
 
 TH2F *hFragp = new TH2F("hFragp","Free spectator protons vs b",50,0.,20.,60,0.,60.);
 hFragp -> SetXTitle("b (fm)");
 hFragp -> SetYTitle("Num. of free spectator protons");
 
 TH2F *hFragn = new TH2F("hFragn","Free spectator neutrons vs b",50,0.,20.,80,0.,80.);
 hFragn -> SetXTitle("b (fm)");
 hFragn -> SetYTitle("Num. of free spectator neutrons");

 TF1 *funb = new  TF1("funb","x",0,20);
 for(Int_t ievent=0; ievent<nev; ievent++){  
   if(ievent%1000==0){printf("Processing event %d\n",ievent);}
   b= Float_t(funb->GetRandom()); 
   spectator(b, &nspectp, &nspectn);
   hsp -> Fill(b,nspectp);
   hsn -> Fill(b,nspectn);
   AliZDCFragment *gallio = new AliZDCFragment(b);
   for(Int_t j=0; j<=99; j++){
      zz[j] =0;
      nn[j] =0;
   }
//
// Generation of fragments
   gallio->GenerateIMF();
   nalpha = gallio->GetNalpha();
   NFrag = gallio->GetFragmentNum();
//
   // Attach neutrons
   ztot = gallio->GetZtot();
   ntot = gallio->GetNtot();
   gallio->AttachNeutrons();
   //
   nspectpfree = (nspectp-ztot-2*nalpha);
   nspectnfree = (nspectn-ntot-2*nalpha);
   
   // Removing deuterons
   ndeu = (Int_t) (nspectnfree*gallio->DeuteronNumber());
   nspectpfree -= ndeu;
   nspectnfree -= ndeu;
   //
   hFragp -> Fill(b, nspectpfree);
   hFragn -> Fill(b, nspectnfree);
//
// Print
   if(debug == 1){
     printf("\n	b = %f",b);
     printf("	#spect p = %d, #spect n = %d\n",nspectp,nspectn); 
     printf("	#spect p free = %d, #spect n free = %d\n",nspectpfree,nspectnfree); 
     printf("	#fragments = %f ",NFrag);
     /*for(Int_t i=0; i<NFrag; i++){
   	printf("\n	ZZ[%d] = %d, NN[%d] = %d\n",i,zz[i],i,nn[i]);
     }*/
     printf("	NAlpha = %d, Ztot = %d, Ntot = %d\n\n",nalpha,ztot,ntot);
   }
   delete gallio;
  } //Event loop
  
  TProfile *profsp = hsp->ProfileX("profsp",-1,-1,"s");
  TProfile *profsn = hsn->ProfileX("profsn",-1,-1,"s");
  TProfile *profFragp = hFragp->ProfileX("profFragp",-1,-1,"s");
  TProfile *profFragn = hFragn->ProfileX("profFragn",-1,-1,"s");

   
  //***********************************************************
  // #### ROOT initialization
  gROOT->Reset();
  gStyle->SetCanvasColor(10);
  gStyle->SetFrameFillColor(10);
  gStyle->SetPalette(1);
  gStyle->SetOptStat(0);
  //
  TCanvas *c1 = new TCanvas("c1","Fragmentation I",0,0,700,700);
  c1->cd();
  c1->SetFillColor(kAzure+7);
  c1->Divide(2,2);
  c1->cd(1);
  hsp -> Draw("colz");   
  c1->cd(2);
  hsn -> Draw("colz");    
  c1->cd(3);
  hFragp -> Draw("colz");   
  c1->cd(4);
  hFragn -> Draw("colz"); 
  //
//.........这里部分代码省略.........