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

void FitterUtils::initiateParams(int nGenSignalZeroGamma, int nGenSignalOneGamma, int nGenSignalTwoGamma, RooRealVar const& expoConstGen, RooRealVar& nSignal, RooRealVar& nPartReco, 
      RooRealVar& nComb, RooRealVar& fracZero, RooRealVar& fracOne, RooRealVar& expoConst, RooRealVar&  nJpsiLeak, bool constPartReco, RooRealVar const& fracPartRecoSigma)
{
   TRandom rand;
   rand.SetSeed();

   int nGenSignal = nGenSignalZeroGamma + nGenSignalOneGamma + nGenSignalTwoGamma;

   double nGenSignal2;
   double nGenPartReco2;
   if(!constPartReco)
   {
      nGenSignal2 = rand.Uniform(nGenSignal-5*sqrt(nGenSignal), nGenSignal+5*sqrt(nGenSignal));
      nGenPartReco2 = rand.Uniform(nGenPartReco-5*sqrt(nGenPartReco), nGenPartReco+5*sqrt(nGenPartReco));
   }
   if(constPartReco)
   { 
      double nGenSigPartReco( nGenSignal+nGenPartReco );
      double nGenSigPartReco2( rand.Uniform( nGenSigPartReco-5*sqrt(nGenSigPartReco), nGenSigPartReco+5*sqrt(nGenSigPartReco) ) );
      double fracPartReco1( nGenPartReco/(1.*nGenSignal));
      double fracPartReco2( rand.Uniform(fracPartReco1-5*fracPartRecoSigma.getVal(), fracPartReco1+5*fracPartRecoSigma.getVal()) ); 

      nGenPartReco2 = fracPartReco2*nGenSigPartReco2 / (1+fracPartReco2); 
      nGenSignal2 = nGenSigPartReco2 / (1+fracPartReco2); 
   }
   double nGenComb2 = rand.Uniform(nGenComb-5*sqrt(nGenComb), nGenComb+5*sqrt(nGenComb));
   double nGenJpsiLeak2 = rand.Uniform(nGenJpsiLeak-5*sqrt(nGenJpsiLeak), nGenJpsiLeak+5*sqrt(nGenJpsiLeak));


   nSignal.setVal(nGenSignal2);
   nSignal.setRange(TMath::Max(0.,nGenSignal2-10.*sqrt(nGenSignal)) , nGenSignal2+10*sqrt(nGenSignal));

   nPartReco.setVal(nGenPartReco2);
   nPartReco.setRange(TMath::Max(0.,nGenPartReco2-10.*sqrt(nGenPartReco)), nGenPartReco2+10*sqrt(nGenPartReco));


   nComb.setVal(nGenComb2);
   nComb.setRange(TMath::Max(0.,nGenComb2-10.*sqrt(nGenComb)), nGenComb2+10*sqrt(nGenComb));

   nJpsiLeak.setVal(nGenJpsiLeak2);
   nJpsiLeak.setRange(TMath::Max(0., nGenJpsiLeak2-10*sqrt(nGenJpsiLeak)), nGenJpsiLeak2+10*sqrt(nGenJpsiLeak));

   double fracGenZero(nGenSignalZeroGamma/(1.*nGenSignal));
   double fracGenOne(nGenSignalOneGamma/(1.*nGenSignal));

   fracZero.setVal(rand.Gaus(fracGenZero, sqrt(nGenSignalZeroGamma)/(1.*nGenSignal))) ;
   fracZero.setRange(0., 1.);
   fracOne.setVal(rand.Gaus(fracGenOne, sqrt(nGenSignalOneGamma)/(1.*nGenSignal))) ;
   fracOne.setRange(0., 1.);

   expoConst.setVal(rand.Uniform( expoConstGen.getVal() - 5*expoConstGen.getError(), expoConstGen.getVal() + 5*expoConstGen.getError() ) );
   expoConst.setRange( expoConstGen.getVal() - 10*expoConstGen.getError(), expoConstGen.getVal() + 10*expoConstGen.getError() );
}

void makeDecalibCDB(Int_t firstRun, Int_t lastRun, Float_t decalib = 0.065)
{
  //Generates a random decalibration factors O(1)
  //to be applied in the anchor run simulations with raw:// .
  //If decalib<0, no decalibration generated, all factors=1.
  
  //Run range is [firstRun,lastRun] and gaussian sigma = decalib.
  //Author: Boris Polishchuk.
  
  AliCDBManager::Instance()->SetDefaultStorage("raw://");
  AliCDBManager::Instance()->SetRun(firstRun);

  TString emcPath("PHOS/Calib/EmcGainPedestals");
  AliCDBEntry* entryEmc = AliCDBManager::Instance()->Get(emcPath.Data(),-1);

  AliPHOSEmcCalibData* clb=0;

  if(entryEmc) clb = (AliPHOSEmcCalibData*)entryEmc->GetObject();
  else { printf("CDB entry not found. Exit.\n"); return; }
  
  if(!clb) { printf("Calibration parameters for PHOS EMC not found.\n"); return; }
  
  printf("\t\tEMC calibration object found: FirstRun=%d LastRun=%d Version=%d SubVersion=%d\n",
         entryEmc->GetId().GetFirstRun(), entryEmc->GetId().GetLastRun(),
         entryEmc->GetId().GetVersion(),entryEmc->GetId().GetSubVersion());
  
  
  TRandom rn;
  rn.SetSeed(0); //the seed is set to the current  machine clock

  Float_t adcChannelEmc;

  for(Int_t module=1; module<6; module++) {
    for(Int_t column=1; column<57; column++) {
      for(Int_t row=1; row<65; row++) {

	if(decalib<0.) adcChannelEmc = 1.;
	else
	  adcChannelEmc =rn.Gaus(1.,decalib);

        clb->SetADCchannelEmcDecalib(module,column,row,adcChannelEmc);
      }
    }
  }

  AliCDBManager::Instance()->SetDefaultStorage("local://./");
  AliCDBStorage* storage = AliCDBManager::Instance()->GetDefaultStorage();
  
  AliCDBMetaData *md = new AliCDBMetaData();
  AliCDBId id(emcPath.Data(),firstRun,lastRun);
  storage->Put(clb,id, md);

}

void JEC_fit_Uncertainty(int N)
{
  TF1 *func[1000];
  TFile *inf = new TFile("L3Graphs_test_Icone5.root");
  TGraphErrors *g = (TGraphErrors*)inf->Get("Correction_vs_CaloPt");
  TGraphErrors *vg[1000];
  int i,k;
  double x[20],y[20],ex[20],ey[20];
  double vx[20],vy[20],vex[20],vey[20];
  for(i=0;i<g->GetN();i++)
    {
      g->GetPoint(i,x[i],y[i]);
      ex[i]=g->GetErrorX(i);
      ey[i]=g->GetErrorY(i); 
    }  
  TRandom *rnd = new TRandom();
  rnd->SetSeed(0);
  for(k=0;k<N;k++)
    {
      for(i=0;i<g->GetN();i++)
        {	
          vx[i] = rnd->Gaus(x[i],ex[i]);
          //vx[i] = x[i];
          vy[i] = rnd->Gaus(y[i],ey[i]);
          vex[i] = ex[i];
          vey[i] = ey[i];
        }
      vg[k] = new TGraphErrors(g->GetN(),vx,vy,vex,vey);
      func[k] = new TF1("func","[0]+[1]/(pow(log10(x),[2])+[3])",1,2000);
      func[k]->SetParameters(1,3,6,5);
      vg[k]->Fit(func[k],"RQ");     	
    }
  
  TCanvas *c = new TCanvas("c","c");
  gPad->SetLogx();
  g->SetMarkerStyle(20);
  g->SetMaximum(3.5);
  g->Draw("AP");
  for(k=0;k<N;k++)
    {
      func[k]->SetLineColor(5);
      func[k]->SetLineWidth(1);
      cout<<func[k]->GetChisquare()<<endl;
      vg[k]->SetMarkerColor(2);
      vg[k]->SetLineColor(2);
      vg[k]->SetMarkerStyle(21);
      //if (func[k]->GetChisquare()<0.1)
        //vg[k]->Draw("sameP");
      func[k]->Draw("same");  
    }  	 
}  

double error(double a, double b) {

 TRandom rndm;
 rndm.SetSeed(12345);
 
 if (a<=0.000000001) return 0.;

 double eff = a/b;
 int n0 = rndm.Poisson((a+b));
 // int n1 = rndm.Binomial(n0,eff);

 double err=sqrt(eff*(1-eff)/n0);

 return err;
}

void fmpmt(void){

  TTree *t = new TTree("tree","fine mesh");
  double adc;
  t->Branch("adc",&adc,"adc/D");

  double Mult(double seed){
    do{
      double _fac=2.79360;
      TRandom _rand;
      _rand.SetSeed(0);
      double _judge = _rand.Uniform(0,1);
      double _x = _rand.Uniform(0.001,10);
    } while (TMath::Gaus(_x, _fac, _fac/20) < _judge);
    return seed*_x;
  }


  double prob = 0.7;
  double npe = 0.1;
  double judge, tmp;
  TRandom rand;
  rand.SetSeed(0);

  for(int i=0; i<4000; ++i){
    if(i%100==0){
      cout<<"=====  "<<i<<"  ====="<<endl;
    }
    do{
      tmp = rand.Uniform(0,4);
      judge = rand.Uniform(0,1);
    } while (TMath::Poisson(tmp, npe) < judge);
    for(int j=0; j<19; ++j){
      if(j==0){
        judge = rand.Uniform(0,1);
        if(judge>prob){
          continue;
        }
      }
      tmp=Mult(tmp);
    }
    adc=tmp;
    t->Fill();
  }
  t->Draw("adc");


}

Int_t foam_kanwa(){
  cout<<"--- kanwa started ---"<<endl;
  TH2D  *hst_xy = new TH2D("hst_xy" ,  "x-y plot", 50,0,1.0, 50,0,1.0);
  Double_t *MCvect =new Double_t[2]; // 2-dim vector generated in the MC run
  //
  TRandom     *PseRan   = new TRandom3();  // Create random number generator
  PseRan->SetSeed(4357);
  TFoam   *FoamX    = new TFoam("FoamX");   // Create Simulator
  FoamX->SetkDim(2);         // No. of dimensions, obligatory!
  FoamX->SetnCells(500);     // Optionally No. of cells, default=2000
  FoamX->SetRhoInt(Camel2);  // Set 2-dim distribution, included below
  FoamX->SetPseRan(PseRan);  // Set random number generator
  FoamX->Initialize();       // Initialize simulator, may take time...
  //
  // visualising generated distribution
  TCanvas *cKanwa = new TCanvas("cKanwa","Canvas for plotting",600,600);
  cKanwa->cd();
  // From now on FoamX is ready to generate events
  int nshow=5000;
  for(long loop=0; loop<100000; loop++){
    FoamX->MakeEvent();            // generate MC event
    FoamX->GetMCvect( MCvect);     // get generated vector (x,y)
    Double_t x=MCvect[0];
    Double_t y=MCvect[1];
    if(loop<10) cout<<"(x,y) =  ( "<< x <<", "<< y <<" )"<<endl;
    hst_xy->Fill(x,y);
    // live plot
    if(loop == nshow){
      nshow += 5000;
      hst_xy->Draw("lego2");
      cKanwa->Update();
    }
  }// loop
  //
  hst_xy->Draw("lego2");  // final plot
  cKanwa->Update();
  //
  Double_t MCresult, MCerror;
  FoamX->GetIntegMC( MCresult, MCerror);  // get MC integral, should be one
  cout << " MCresult= " << MCresult << " +- " << MCerror <<endl;
  cout<<"--- kanwa ended ---"<<endl;

  return 0;
}//kanwa

void ThrowTheDices() {

  TH1S* dices = new TH1S("dices","Throw with 3 dices;Pips;Count rate", 20,0.5,20.5);

  TRandom* die = new TRandom();
  die->SetSeed(0);

  Float_t* value = new Float_t[3];

  for (Int_t i=0; i<10000;++i) {
    die->RndmArray(3,value);
    Int_t sum = static_cast<Int_t>((value[0]*6) + 1) +
                static_cast<Int_t>((value[1]*6) + 1) +
                static_cast<Int_t>((value[2]*6) + 1);
    dices->Fill(sum);
  }

  dices->Draw();

}

void FitterUtilsSimultaneousExpOfPolyTimesX::initiateParams(int nGenSignalZeroGamma, int nGenSignalOneGamma, int nGenSignalTwoGamma,
      RooRealVar& nKemu, RooRealVar& nSignal, RooRealVar& nPartReco,
      RooRealVar& nComb, RooRealVar& fracZero, RooRealVar& fracOne,
      RooRealVar&  nJpsiLeak, bool constPartReco, RooRealVar const& fracPartRecoSigma,
      RooRealVar& l1Kee, RooRealVar& l2Kee, RooRealVar& l3Kee, RooRealVar& l4Kee, RooRealVar& l5Kee,
      RooRealVar& l1Kemu, RooRealVar& l2Kemu, RooRealVar& l3Kemu, RooRealVar& l4Kemu, RooRealVar& l5Kemu,
      RooRealVar const& l1KeeGen, RooRealVar const& l2KeeGen, RooRealVar const& l3KeeGen, RooRealVar const& l4KeeGen, RooRealVar const& l5KeeGen 

      )
{
  FitterUtilsExpOfPolyTimesX::initiateParams(nGenSignalZeroGamma, nGenSignalOneGamma, nGenSignalTwoGamma,
           nSignal, nPartReco, nComb, fracZero, fracOne, nJpsiLeak, constPartReco, fracPartRecoSigma,
           l1Kee, l2Kee, l3Kee, l4Kee, l5Kee,
           l1KeeGen, l2KeeGen, l3KeeGen, l4KeeGen, l5KeeGen ); 



  TRandom rand;
  rand.SetSeed();

  nKemu.setVal(rand.Uniform(nGenKemu-5*sqrt(nGenKemu), nGenKemu+5*sqrt(nGenKemu)));
  nKemu.setRange(nGenKemu-10*sqrt(nGenKemu), nGenKemu+10*sqrt(nGenKemu));

  l1Kemu.setVal(rand.Uniform( l1KeeGen.getVal() - 5*l1KeeGen.getError(), l1KeeGen.getVal() + 5*l1KeeGen.getError() ) );
  l1Kemu.setRange( l1KeeGen.getVal() - 10*l1KeeGen.getError(), l1KeeGen.getVal() + 10*l1KeeGen.getError() );

  l2Kemu.setVal(rand.Uniform( l2KeeGen.getVal() - 5*l2KeeGen.getError(), l2KeeGen.getVal() + 5*l2KeeGen.getError() ) );
  l2Kemu.setRange( l2KeeGen.getVal() - 10*l2KeeGen.getError(), l2KeeGen.getVal() + 10*l2KeeGen.getError() );

  l3Kemu.setVal(rand.Uniform( l3KeeGen.getVal() - 5*l3KeeGen.getError(), l3KeeGen.getVal() + 5*l3KeeGen.getError() ) );
  l3Kemu.setRange( l3KeeGen.getVal() - 10*l3KeeGen.getError(), l3KeeGen.getVal() + 10*l3KeeGen.getError() );

  l4Kemu.setVal(rand.Uniform( l4KeeGen.getVal() - 5*l4KeeGen.getError(), l4KeeGen.getVal() + 5*l4KeeGen.getError() ) );
  l4Kemu.setRange( l4KeeGen.getVal() - 10*l4KeeGen.getError(), l4KeeGen.getVal() + 10*l4KeeGen.getError() );

  l5Kemu.setVal(rand.Uniform( l5KeeGen.getVal() - 5*l5KeeGen.getError(), l5KeeGen.getVal() + 5*l5KeeGen.getError() ) );
  l5Kemu.setRange( l5KeeGen.getVal() - 10*l5KeeGen.getError(), l5KeeGen.getVal() + 10*l5KeeGen.getError() );

}

void FitBias(TString CAT,TString CUT,float SIG,float BKG,int NTOYS)
{
  gROOT->ForceStyle();
  
  RooMsgService::instance().setSilentMode(kTRUE);
  RooMsgService::instance().setStreamStatus(0,kFALSE);
  RooMsgService::instance().setStreamStatus(1,kFALSE);
  
  // -----------------------------------------
  TFile *fTemplates = TFile::Open("templates_"+CUT+"_"+CAT+"_workspace.root");
  RooWorkspace *wTemplates = (RooWorkspace*)fTemplates->Get("w");
  RooRealVar *x            = (RooRealVar*)wTemplates->var("mTop");
  RooAbsPdf *pdf_signal    = (RooAbsPdf*)wTemplates->pdf("ttbar_pdf_Nominal");
  RooAbsPdf *pdf_bkg       = (RooAbsPdf*)wTemplates->pdf("qcdCor_pdf"); 
  TRandom *rnd = new TRandom();
  rnd->SetSeed(0);
  x->setBins(250);   
  RooPlot *frame;

  TFile *outf;

  if (NTOYS > 1) { 
    outf = TFile::Open("FitBiasToys_"+CUT+"_"+CAT+".root","RECREATE");
  }

  float nSigInj,nBkgInj,nSigFit,nBkgFit,eSigFit,eBkgFit,nll;

  TTree *tr = new TTree("toys","toys");
  
  tr->Branch("nSigInj",&nSigInj,"nSigInj/F");
  tr->Branch("nSigFit",&nSigFit,"nSigFit/F");
  tr->Branch("nBkgInj",&nBkgInj,"nBkgInj/F");
  tr->Branch("nBkgFit",&nBkgFit,"nBkgFit/F");
  tr->Branch("eSigFit",&eSigFit,"eSigFit/F");
  tr->Branch("eBkgFit",&eBkgFit,"eBkgFit/F");
  tr->Branch("nll"    ,&nll    ,"nll/F");

  for(int itoy=0;itoy<NTOYS;itoy++) {
    // generate pseudodataset
    nSigInj = rnd->Poisson(SIG);
    nBkgInj = rnd->Poisson(BKG);
    RooRealVar *nSig = new RooRealVar("nSig","nSig",nSigInj);
    RooRealVar *nBkg = new RooRealVar("nBkg","nBkg",nBkgInj);
    RooAddPdf *model = new RooAddPdf("model","model",RooArgList(*pdf_signal,*pdf_bkg),RooArgList(*nSig,*nBkg)); 
    RooDataSet *data = model->generate(*x,nSigInj+nBkgInj);
    
    RooDataHist *roohist = new RooDataHist("roohist","roohist",RooArgList(*x),*data);
    // build fit model
    RooRealVar *nFitSig = new RooRealVar("nFitSig","nFitSig",SIG,0,10*SIG);
    RooRealVar *nFitBkg = new RooRealVar("nFitBkg","nFitBkg",BKG,0,10*BKG);
    RooAddPdf *modelFit = new RooAddPdf("modelFit","modelFit",RooArgList(*pdf_signal,*pdf_bkg),RooArgList(*nFitSig,*nFitBkg)); 
    // fit the pseudo dataset
    RooFitResult *res = modelFit->fitTo(*roohist,RooFit::Save(),RooFit::Extended(kTRUE));
    //res->Print();
    nSigFit = nFitSig->getVal();
    nBkgFit = nFitBkg->getVal();
    eSigFit = nFitSig->getError();
    eBkgFit = nFitBkg->getError();
    nll     = res->minNll();
    tr->Fill();
    if (itoy % 100 == 0) {
      cout<<"Toy #"<<itoy<<": injected = "<<nSigInj<<", fitted = "<<nSigFit<<", error = "<<eSigFit<<endl;
    }
    if (NTOYS == 1) {
      frame = x->frame();
      roohist->plotOn(frame); 
      model->plotOn(frame);
    }
  }
  if (NTOYS == 1) {
    TCanvas *can = new TCanvas("Toy","Toy",900,600);
    frame->Draw();
  }  
  else {
    outf->cd();
    tr->Write();
    outf->Close();
    fTemplates->Close();
  }  
}

//.........这里部分代码省略.........
  FitResults fr_pedSubtracted_0 = fitSingleHisto( h1_cef3_pedSubtracted_0, 0., 0., lowerFitBoundary[0]-pedestals[0].first, 185.-pedestals[0].first );
  FitResults fr_pedSubtracted_1 = fitSingleHisto( h1_cef3_pedSubtracted_1, 0., 0., lowerFitBoundary[1]-pedestals[1].first, 190.-pedestals[1].first );
  FitResults fr_pedSubtracted_2 = fitSingleHisto( h1_cef3_pedSubtracted_2, 0., 0., lowerFitBoundary[2]-pedestals[2].first, 190.-pedestals[2].first );
  FitResults fr_pedSubtracted_3 = fitSingleHisto( h1_cef3_pedSubtracted_3, 0., 0., lowerFitBoundary[3]-pedestals[3].first, 198.-pedestals[3].first );


  std::vector<float> correctionFactors = intercalibrateFibers(fr_0,fr_1,fr_2,fr_3,true);
  std::vector<float> correctionFactors_pedSubtracted = intercalibrateFibers(fr_pedSubtracted_0,fr_pedSubtracted_1,fr_pedSubtracted_2,fr_pedSubtracted_3,true);
  std::vector<float> correctionFactors_pedSubtracted_muQ = intercalibrateFibers(fr_pedSubtracted_0,fr_pedSubtracted_1,fr_pedSubtracted_2,fr_pedSubtracted_3,false);

  std::vector<float> mu;
  mu.push_back(fr_0.mu);
  mu.push_back(fr_1.mu);
  mu.push_back(fr_2.mu);
  mu.push_back(fr_3.mu);
  float mumean=0.;


  for(int i=0;i<4;i++){
    mumean+=mu[i];
  }

  mumean=mumean/4.;

  std::vector<float> correctionFactors_pedSubtracted_muMean= intercalibrateFibers(fr_pedSubtracted_0,fr_pedSubtracted_1,fr_pedSubtracted_2,fr_pedSubtracted_3,false,mumean);

  for(int i=0;i<4;i++){
    std::cout<<"correctionFactors "<<correctionFactors[i]<<std::endl;
    std::cout<<"correctionFactors pedSub "<<correctionFactors_pedSubtracted[i]<<std::endl;
    std::cout<<"correctionFactors pedSub muQ "<<correctionFactors_pedSubtracted_muQ[i]<<std::endl;
    std::cout<<"correctionFactors pedSub muMean "<<correctionFactors_pedSubtracted_muMean[i]<<std::endl;
  }
  TRandom a;
  a.SetSeed(100);

  for( unsigned iEntry=0; iEntry<nentries; ++iEntry ) {
    
    tree->GetEntry(iEntry);
    
    if( iEntry % 5000 == 0 ) std::cout << "Entry: " << iEntry << " / " << nentries << std::endl;
    
    for( unsigned i=0; i<40; ++i ) {
      
      int board  = adcBoard[i];
      int channel= adcChannel[i];
      
      
      float cef3_corr=0;

      
      if( board==CEF3_ADC_BOARD ) {
	if( channel==(CEF3_ADC_START_CHANNEL  ) ){
	  h1_cef3_corr_0->Fill(adcData[i]*correctionFactors[0]);
	  if(adcData[i]>(pedestals[0].first + nSigma*pedestals[0].second)){
	    h1_cef3_pedSubtracted_corr_0->Fill((adcData[i]-pedestals[0].first)*correctionFactors_pedSubtracted[0]);
	    h1_cef3_pedSubtracted_corr_rebin_0->Fill((adcData[i]-pedestals[0].first)*correctionFactors_pedSubtracted[0]);
	    h1_cef3_pedSubtracted_corr_muQ_0->Fill((adcData[i]-pedestals[0].first)*correctionFactors_pedSubtracted_muQ[0]);
	    h1_cef3_pedSubtracted_corr_muMean_0->Fill((adcData[i]-pedestals[0].first)*correctionFactors_pedSubtracted_muMean[0]);
	    h1_cef3_pedSubtracted_corr_sum->Fill(a.Uniform(adcData[i]-pedestals[0].first-0.5,adcData[i]-pedestals[0].first+0.5)*correctionFactors_pedSubtracted[0]);
	  }
	  cef3_corr+=adcData[i]*correctionFactors[0];
	}
        else if( channel==(CEF3_ADC_START_CHANNEL+1) ){
	  h1_cef3_corr_1->Fill(adcData[i]*correctionFactors[1]);
	  if(adcData[i]>(pedestals[1].first + nSigma*pedestals[1].second)) {
	    h1_cef3_pedSubtracted_corr_1->Fill((adcData[i]-pedestals[1].first)*correctionFactors_pedSubtracted[1]);

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

  //std::cout << "seedMin..seedMax=" << seedMin << ".." << seedMax << "\n";

  //return retCodeOk;

  //--------------------------------------------------------------------------------------------------------------
  // Main analysis code
  //==============================================================================================================


  std::cout << mainpart;

  TRandom random;
  std::vector<ElectronEnergyScale*> escaleV;
  std::vector<EventWeight_t*> specEWeightsV;
  std::vector<double> specReweightsV;
  std::vector<EventSelector_t*> evtSelectorV;
  std::vector<TH2D*> specTH2DWeightV; // used for ESCALE_RESIDUAL

  double specWeight=1.;
  int useSpecWeight=0;
  if (systMode==DYTools::FSR_5plus) { specWeight=1.05; useSpecWeight=1; }
  else if (systMode==DYTools::FSR_5minus) { specWeight=0.95; useSpecWeight=1; }
  else if (systMode==DYTools::FSR_RND_STUDY) useSpecWeight=1;

  // check random seed. Special weights use their own,
  // built-in dependencies on seed
  {
    int startSeed=-1;
    if (systMode==DYTools::SYST_RND) {
      std::cout << "setting startSeed=" << globalSeed << "\n";
      startSeed= globalSeed;
    }
    random.SetSeed(startSeed);
    gRandom->SetSeed(startSeed);
  }

  // The random seeds are needed only if we are running this script in systematics mode

  if (systMode==DYTools::FSR_STUDY) {
    specReweightsV.reserve(seedDiff);
    specEWeightsV.reserve(seedDiff);
    for (int i=0; i<seedDiff; ++i) {
      double specW= 1 + 0.05*(i-1);
      specReweightsV.push_back(specW);
      specEWeightsV.push_back(new EventWeight_t(evWeight));
    }
  }
  else if (systMode==DYTools::PU_STUDY) {
    if (inpMgr.puReweightFlag()==0) {
      std::cout << "systMode=PU_STUDY needs puReweightFlag=1 in the input file\n";
      return retCodeError;
    }
    specEWeightsV.reserve(2);
    for (int i=0; i<2; ++i) {
      DYTools::TSystematicsStudy_t study=(i==0) ? DYTools::PILEUP_5minus : DYTools::PILEUP_5plus;
      EventWeight_t *ew=new EventWeight_t();
      if (!ew->init(inpMgr.puReweightFlag(),inpMgr.fewzFlag(),study,rndStudyStr)) {
	std::cout << "in plotUnfoldingMatrix.C\n";
	return retCodeError;
      }
      specEWeightsV.push_back(ew);
    }
  }
  else if (systMode==DYTools::SYST_RND) {
    // nothing special about weights

bool PSEC4Sim::Execute(){
  
  PSECElement *tmp=new PSECElement(6);
  
  tmp->starttime=0;
  tmp->timestep=acq_rate_usecs;
  tmp->numpoints=256;
  tmp->numreadouts=num_readouts;

  m_data->vars.Get("crateid",tmp->crateid);
  m_data->vars.Get("cardid",tmp->cardid);
  tmp->cardid++;
  if(tmp->cardid>7){
    tmp->crateid++;
    tmp->cardid=0;
    m_data->vars.Set("crateid",tmp->crateid);
  }
  m_data->vars.Set("cardid",tmp->cardid);
  std::cout<<"cardid "<<tmp->cardid<<" crateid "<<tmp->crateid<<std::endl;
  int tmppmtid;
  m_data->vars.Get("pmtid",tmppmtid);
  tmppmtid+=1;
  m_data->vars.Set("pmtid",tmppmtid+5);

  double pi=3.14159265351;
  double angle =(2*pi)/(6*8);
  double r=5.5;
  double dist=22.0/72.0;
  

  for (int i=0;i<6;i++){
    tmp->pmtids[i]=i+tmppmtid;
    std::vector<float> tmploc;
    tmploc.push_back(cos((((tmp->cardid)*6)+(i+1))*angle)*r);
     std::cout<<"i= "<<i<<" angle multiplier "<<((tmp->cardid)*6)+(i+1)<<std::endl;

    tmploc.push_back(sin((((tmp->cardid)*6)+(i+1))*angle)*r);
    tmploc.push_back((dist/2.0)+(tmp->crateid*dist));
    tmploc.push_back((((tmp->cardid)*6)+(i+1))*angle);
    tmp->pmtlocations.push_back(tmploc);
  }
  
  
  TRandom *eventGenerator = new TRandom();
  
time_t seconds;

// seconds = time (NULL);
 uint seed= ((dist/2.0)+(tmp->crateid*dist)+((tmp->cardid)*6))*100;
  //std::cout<<"seed = "<<((uint) seconds)<<std::endl;
  eventGenerator->SetSeed(seed);
  std::vector<std::vector<std::vector<float> > >tmpf3;
  std::vector< std::vector<std::vector<int> > >tmpi3;


  for(int k=0;k<num_readouts;k++){

  std::vector<std::vector<float> >tmpf2;
  std::vector<std::vector<int> >tmpi2;
  
  for(int j=0;j<6;j++){
    std::vector<float> tmpf;
    std::vector<int> tmpi;
    
    for (int i=0;i<256;i++){
      tmpf.push_back(eventGenerator->Gaus(0.00,0.002));
      tmpi.push_back(floor(eventGenerator->Gaus(2100.0,50.0)));
    }
    tmpf2.push_back(tmpf);
    tmpi2.push_back(tmpi);
  }

  tmpf3.push_back(tmpf2);
  tmpi3.push_back(tmpi2);
  }
  
  tmp->caldata=tmpf3;
  tmp->rawdata=tmpi3;
  
  m_data->elements.push_back(tmp);

  return true;

}

//.........这里部分代码省略.........
	    return;
	  }
	}
	std::cout << "energy scale corrections: " << escale.calibrationSetName() << "\n";
      }
      string fname;
      Int_t color, linesty;
      stringstream ss(line);
      Double_t xsec;
      ss >> fname >> xsec >> color >> linesty;
      string label = line.substr(line.find('@')+1);
      fnamev.push_back(fname);
      labelv.push_back(label);
      colorv.push_back(color);
      linev.push_back(linesty);
      xsecv.push_back(xsec);
      lumiv.push_back(0);
    }
  }
  ifs.close();
  
  const Double_t kGAP_LOW  = 1.4442;
  const Double_t kGAP_HIGH = 1.566;


  //--------------------------------------------------------------------------------------------------------------
  // Main analysis code 
  //==============================================================================================================


  TRandom random;
  // The random seeds are needed only if we are running this script in systematics mode
  int seed = randomSeed;
  random.SetSeed(seed);
  gRandom->SetSeed(seed);
  // In the case of systematic studies, generate an array of random offsets
  TVectorD shift(escale._nEtaBins); // this vector is outdated by the new features in the escale obj.class
  shift = 0;
  if(systematicsMode==DYTools::RESOLUTION_STUDY) {
    escale.randomizeSmearingWidth(seed);
    for(int i=0; i<escale._nEtaBins; i++)
      shift[i] = gRandom->Gaus(0,1);
  }

  // prepare tools for ESCALE_RESIDUAL
  TH1F *shapeWeights=NULL;
  if (systematicsMode==DYTools::ESCALE_RESIDUAL) {
    TString shapeFName=TString("../root_files/yields/") + dirTag + TString("/shape_weights.root");
    std::cout << "Obtaining shape_weights.root from <" << shapeFName << ">\n";
    TFile fshape(shapeFName);
    if (!fshape.IsOpen()) {
      std::cout << "failed to open a file <" << shapeFName << ">\n";
      throw 2;
    }
    shapeWeights = (TH1F*)fshape.Get("weights");
    shapeWeights->SetDirectory(0);
    dirTag += TString("_escale_residual");
    std::cout << "changing dirTag to <" << dirTag << ">\n";
  }

  //  
  // Set up histograms
  //
  vector<TH1F*> hZMassv;//, hZMass2v, hZPtv, hZPt2v, hZyv, hZPhiv;  
  
  char hname[100];

int main(int argc, char *argv[])
{
    TRandom rng;
    rng.SetSeed(2014);
    hxx_tree data;

    std::string opt, infile, outroot, outdir;
    koptions options(argc, argv);

    //check for the --help option:
    if ( options.find("--help") ) {
        usage();
    }
    //options.set("--seed=", seed);
    //if (seed > 0) rng.SetSeed(seed);
    double met_smear = 30.0;
    options.set("--met_smear=", met_smear);
    int    num_smear = 1;
    options.set("--num_smear=", num_smear);
    double fake_rate = 0.0;
    options.set("--fake_rate=", fake_rate);
    int mode_8tev = options.find("--8tev");
    int pub_plots = options.find("--pub_plots");



    double lumi = 300;
    double met_cut = 250.0;

    if (mode_8tev) {
        lumi = 20.0;
        met_cut = 250;
        cout << "INFO:  settings for sqrt(s) = 8 TeV\n";
    } else {
        cout << "INFO:  settings for sqrt(s) = 14 TeV\n";
    }
    cout << "INFO:  lumi:     " << lumi << "\n";
    cout << "INFO:  MET cut:  " << met_cut << "\n";

    //check for unrecognized options (beginning with -- or -)
    while(options.unrecognized_option(opt)) {
        cout << "WARNING: unrecognized option:" << opt << "\n";
        usage();
    }

    if (options.size() != 3) {
        usage();
    }

    options.shift_argument(infile);
    options.shift_argument(outroot);
    options.shift_argument(outdir);

    cout << "INFO: reading analysis tree file:    " << infile  << "\n";
    cout << "INFO: writing histogram root file:   " << outroot << "\n";
    cout << "INFO: writing results to directory:  " << outdir  << "\n";
    cout << "INFO: MET smearing amount is " << met_smear << "\n";
    cout << "INFO: MET smearing number is " << num_smear << "\n";
    cout << "INFO: W+jets fake rate is " << fake_rate << "\n";

    auto_write aw;

    cutflow_tool cutflow;
    histogram_manager h0mll(new TH1F("h0mll","",60,60.0,120.0));

    if (pub_plots) {
        h0mll.add_sample(1,  "_z");
        h0mll.add_sample(2,  "_h");
        h0mll.add_sample(3,  "_w");
        h0mll.add_sample(4,  "_top");

        cutflow.add_sample_name(1, "Z,ZZ,ZW");
        cutflow.add_sample_name(2, "Higgs");
        cutflow.add_sample_name(3, "W,WW");
        cutflow.add_sample_name(4, "tt");
    } else {
        h0mll.add_sample(1,  "_zjj");
        h0mll.add_sample(2,  "_zz_zw");
        h0mll.add_sample(3,  "_ww");
        h0mll.add_sample(4,  "_tt");
        h0mll.add_sample(5,  "_hzz");
        h0mll.add_sample(6,  "_zh");
        h0mll.add_sample(7,  "_wh");
        h0mll.add_sample(8,  "_wjjj");

        cutflow.add_sample_name(1, "Zjj");
        cutflow.add_sample_name(2, "ZZ,ZW");
        cutflow.add_sample_name(3, "WW");
        cutflow.add_sample_name(4, "tt");
        cutflow.add_sample_name(5, "HZZ");
        cutflow.add_sample_name(6, "ZH");
        cutflow.add_sample_name(7, "WH");
        cutflow.add_sample_name(8, "W+jets");
    }




    h0mll.add_sample(20,  "_hxx1");
    h0mll.add_sample(21,  "_hxx10");
//.........这里部分代码省略.........