C++ TStopwatch::Start方法代码示例

int main(int argc, char ** argv) {
	TStopwatch comulativeWatch;
	comulativeWatch.Start(true);
	readInputs(argc,argv);
	cout<<"Currrent Subversion Revision: "<<SVN_REV<<endl;
	cout << "starting main loop.." << endl;
	RunListOK = ReadRunList();
	if(!RunListOK)exit(-1);
	TSystem* sys = gSystem;
	std::string currentDir = sys->pwd();
	for (unsigned int i = 0; i < RunParameters.size(); i++) {
		cout << RunParameters[i].getRunNumber();
		if (i+1 < RunParameters.size()) cout << ", ";
	}
	cout << " will be analysed.." << endl;

	if (!RunListOK) return 0;

	/**Start with Analyising, read RunParameteres of the Run and start analysis with that parameters
	 */
	for (unsigned int i = 0; i < RunParameters.size(); i++) {
		TStopwatch runWatch;
		runWatch.Start(true);
		RunParameters[i].Print();
		bool DO_ALIGNMENT = RunParameters[i].doAlignment();
		bool DO_ALIGNMENTANALYSIS = RunParameters[i].doAlignmentAnalysis();
		bool DO_TRANSPARENT_ANALYSIS = RunParameters[i].doTransparentAnalysis();

		time_t rawtime;
		tm *timestamp;
		time (&rawtime);
		timestamp = gmtime(&rawtime);

		ostringstream logfilename;
		logfilename << "analyse_log_" << RunParameters[i].getRunNumber() << "_" << timestamp->tm_year << "-" << timestamp->tm_mon << "-" << timestamp->tm_mday << "." << timestamp->tm_hour << "." << timestamp->tm_min << "." << timestamp->tm_sec << ".log";
		//
		//		FILE *log;
		//		log = freopen(logfilename.str().c_str(), "w", stdout);


		TSettings *settings = 0;
		cout<<"settings"<<endl;
		settings = new TSettings((TRunInfo*)&RunParameters[i]);

		TResults *currentResults =new TResults(settings);
		currentResults->Print();


		TRawEventSaver *eventSaver;
		eventSaver = new TRawEventSaver(settings);
		eventSaver->saveEvents(RunParameters[i].getEvents());
		delete eventSaver;

		//Calculate Pedestal
		sys->cd(currentDir.c_str());
		TPedestalCalculation* pedestalCalculation;
		pedestalCalculation = new TPedestalCalculation(settings);
		pedestalCalculation->calculateSlidingPedestals(RunParameters[i].getEvents());
		delete pedestalCalculation;

		if(RunParameters[i].doPedestalAnalysis()){
			sys->cd(currentDir.c_str());
			TAnalysisOfPedestal *analysisOfPedestal;
			analysisOfPedestal = new TAnalysisOfPedestal(settings);
			analysisOfPedestal->setResults(currentResults);
			analysisOfPedestal->doAnalysis(RunParameters[i].getEvents());
			delete analysisOfPedestal;
		}

		THTMLGenerator *htmlGen = new THTMLGenerator(settings);
		htmlGen->setFileGeneratingPath(settings->getAbsoluteOuputPath(true));
		htmlGen->setMainPath("./");
		htmlGen->setSubdirPath("");
		htmlGen->setFileName("overview.html");
		htmlGen->addSection("Pedestal","<a href=\"./pedestalAnalysis/pedestal.html\">PEDESTAL</a>");
		htmlGen->addSection("Clustering","<a href=\"./clustering/clustering.html\">CLUSTERING</a>");
		htmlGen->addSection("Selection","<a href=\"./selections/selection.html\">SELECTION</a>");
		htmlGen->addSection("Alignment","<a href=\"./alignment/alignment.html\">ALIGNMENT</a>");
		htmlGen->addSection("Landau","<a href=\"./selectionAnalysis/landaus.html\">LANDAU-DISTRIBUTIONS</a>");
		htmlGen->generateHTMLFile();
		delete htmlGen;

		sys->cd(currentDir.c_str());
		TClustering* clustering;
		clustering=new TClustering(settings);//int seedDetSigma=10,int hitDetSigma=7,int seedDiaSigma=5, int hitDiaSigma=3);
		clustering->ClusterEvents(RunParameters[i].getEvents());
		delete clustering;

		if(RunParameters[i].doClusterAnalysis()){
			sys->cd(currentDir.c_str());
			TAnalysisOfClustering* analysisClustering;
			analysisClustering = new TAnalysisOfClustering(settings);
			analysisClustering->setResults(currentResults);
			analysisClustering->doAnalysis(RunParameters[i].getEvents());
			delete analysisClustering;
		}

		sys->cd(currentDir.c_str());
		TSelectionClass* selectionClass = new TSelectionClass(settings);
		selectionClass->SetResults(currentResults);
//.........这里部分代码省略.........

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

  // load framework libraries
  gSystem->Load( "libFWCoreFWLite" );
  AutoLibraryLoader::enable();

  if ( argc < 2 ) {
    std::cout << "Usage : " << argv[0] << " [parameters.py]" << std::endl;
    return 0;
  }

  // Get the python configuration
  PythonProcessDesc builder(argv[1], argc, argv);
  edm::ParameterSet const& shyftParameters = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("shyftAnalysis");
  edm::ParameterSet const& inputs = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("inputs");
  edm::ParameterSet const& outputs = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("outputs");

  // book a set of histograms
  fwlite::TFileService fs = fwlite::TFileService( outputs.getParameter<std::string>("outputName") );
  TFileDirectory theDir = fs.mkdir( "histos" ); 
    
  // This object 'event' is used both to get all information from the
  // event as well as to store histograms, etc.
  fwlite::ChainEvent ev ( inputs.getParameter<std::vector<std::string> > ("fileNames") );

  //cout << "Making event selector" << endl;
  WPlusJetsEventSelector wPlusJets( shyftParameters );
  pat::strbitset ret = wPlusJets.getBitTemplate();
  

  unsigned int nEventsAnalyzed(0);

  // some timing
  TStopwatch timer;
  timer.Start();

  std::cout << "About to loop" << std::endl;
  //loop through each event
  for( ev.toBegin();
       ! ev.atEnd();
       ++ev) {
    ret.set(false);
    //bool passed = wPlusJets(ev, ret);

    ++nEventsAnalyzed;

  } //end event loop

  timer.Stop();

  // print some timing statistics
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
  printf("Analyzed events: %d \n",nEventsAnalyzed);
  printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
  printf("%4.2f events / RealTime second .\n", nEventsAnalyzed/rtime);
  printf("%4.2f events / CpuTime second .\n", nEventsAnalyzed/ctime);


  
  //cout << "Printing" << endl;
  wPlusJets.print(std::cout);
  //cout << "We're done!" << endl;

  
  return 0;
}

void compute_p0(const char* inFileName,
	    const char* wsName = "combined",
	    const char* modelConfigName = "ModelConfig",
	    const char* dataName = "obsData",
	    const char* asimov1DataName = "asimovData_1",
	    const char* conditional1Snapshot = "conditionalGlobs_1",
	    const char* nominalSnapshot = "nominalGlobs",
	    string smass = "130",
	    string folder = "test")
{
  double mass;
  stringstream massStr;
  massStr << smass;
  massStr >> mass;

  double mu_profile_value = 1; // mu value to profile the obs data at wbefore generating the expected
  bool doConditional      = 1; // do conditional expected data
  bool remakeData         = 0; // handle unphysical pdf cases in H->ZZ->4l
  bool doUncap            = 1; // uncap p0
  bool doInj              = 0; // setup the poi for injection study (zero is faster if you're not)
  bool doObs              = 1; // compute median significance
  bool doMedian           = 1; // compute observed significance

  TStopwatch timer;
  timer.Start();

  TFile f(inFileName);
  RooWorkspace* ws = (RooWorkspace*)f.Get(wsName);
  if (!ws)
  {
    cout << "ERROR::Workspace: " << wsName << " doesn't exist!" << endl;
    return;
  }
  ModelConfig* mc = (ModelConfig*)ws->obj(modelConfigName);
  if (!mc)
  {
    cout << "ERROR::ModelConfig: " << modelConfigName << " doesn't exist!" << endl;
    return;
  }
  RooDataSet* data = (RooDataSet*)ws->data(dataName);
  if (!data)
  {
    cout << "ERROR::Dataset: " << dataName << " doesn't exist!" << endl;
    return;
  }

  mc->GetNuisanceParameters()->Print("v");

  ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2");
  ROOT::Math::MinimizerOptions::SetDefaultStrategy(0);
  ROOT::Math::MinimizerOptions::SetDefaultPrintLevel(1);
  cout << "Setting max function calls" << endl;

  ws->loadSnapshot("conditionalNuis_0");
  RooArgSet nuis(*mc->GetNuisanceParameters());

  RooRealVar* mu = (RooRealVar*)mc->GetParametersOfInterest()->first();

  RooAbsPdf* pdf = mc->GetPdf();

  string condSnapshot(conditional1Snapshot);
  RooArgSet nuis_tmp2 = *mc->GetNuisanceParameters();
  RooNLLVar* obs_nll = doObs ? (RooNLLVar*)pdf->createNLL(*data, Constrain(nuis_tmp2)) : NULL;

  RooDataSet* asimovData1 = (RooDataSet*)ws->data(asimov1DataName);
  RooRealVar* emb = (RooRealVar*)mc->GetNuisanceParameters()->find("ATLAS_EMB");
  if (!asimovData1 || (string(inFileName).find("ic10") != string::npos && emb))
  {
    if (emb) emb->setVal(0.7);
    cout << "Asimov data doesn't exist! Please, allow me to build one for you..." << endl;
    string mu_str, mu_prof_str;
    asimovData1 = makeAsimovData(mc, doConditional, ws, obs_nll, 1, &mu_str, &mu_prof_str, mu_profile_value, true);
    condSnapshot="conditionalGlobs"+mu_prof_str;
  }
  
  if (!doUncap) mu->setRange(0, 40);
  else mu->setRange(-40, 40);

  RooAbsPdf* pdf = mc->GetPdf();

  RooArgSet nuis_tmp1 = *mc->GetNuisanceParameters();
  RooNLLVar* asimov_nll = (RooNLLVar*)pdf->createNLL(*asimovData1, Constrain(nuis_tmp1));

  //do asimov
  mu->setVal(1);
  mu->setConstant(0);
  if (!doInj) mu->setConstant(1);

  int status,sign;
  double med_sig=0,obs_sig=0,asimov_q0=0,obs_q0=0;

  if (doMedian)
  {
    ws->loadSnapshot(condSnapshot.c_str());
    if (doInj) ws->loadSnapshot("conditionalNuis_inj");
    else ws->loadSnapshot("conditionalNuis_1");
    mc->GetGlobalObservables()->Print("v");
    mu->setVal(0);
    mu->setConstant(1);
    status = minimize(asimov_nll, ws);
//.........这里部分代码省略.........

void DrawFakeEtaBins_reduced(std::string kSpecies="PbPb"){
  timer.Start();

  bool printDebug=false;
   TDatime date;

  std::string kAlgName="ak";
  std::string kSpeciesLong = "pp (PYTHIA)";
  if( kSpecies == "PbPb" ) { kAlgName="akPu";  kSpeciesLong = "PbPb (PYTHIA+HYDJET)"; }
  std::string kjetType="PF";
  
  LoadStyle();
  TH1::SetDefaultSumw2();
  TH2::SetDefaultSumw2();
  TH2F * hBlankRatio = new TH2F("hBlankRatio","",20,-2.1,2.1,50,0,1.0);
  TH2F * hBlankPhiRatio = new TH2F("hBlankPhiRatio","",20,-pi,pi,50,0,1.0);
  TH2F * hBlankPt= new TH2F("hBlankPt","",50,20,350,50,0,1.0);
  TH2F * hBlankPtEff= new TH2F("hBlankPtEff","",50,20,350,50,0.6,1.1);


    TH1F *hPtAll [2][ncen], *hPtEff [2][ncen];
    TH1F *hPtFakeAll[ncen], *hPtFake[2][ncen];
    TH1F *hPtFakeRatio[2][ncen];
    TH1F *hPtAll_etabin[2][ncen][neta], *hPtEff_etabin[2][ncen][neta];    
    TH1F *hPtFakeAll_etabin[ncen][neta], *hPtFake_etabin[2][ncen][neta];
    TH1F *hPtFakeRatio_etabin[2][ncen][neta];   
    TH1F *hPtEffRatio_etabin[2][ncen][neta];     
    TH1F *hPtEffRatio[2][ncen];     
    TH1F *hEtaFakeAll_20[ncen],*hEtaFakeAll_30[ncen], *hEtaFakeAll_40[ncen],*hEtaFakeAll_45[ncen], *hEtaFakeAll_50[ncen], *hEtaFakeAll_60[ncen];
    TH1F *hPhiFakeAll_20[ncen],*hPhiFakeAll_30[ncen], *hPhiFakeAll_40[ncen],*hPhiFakeAll_45[ncen], *hPhiFakeAll_50[ncen], *hPhiFakeAll_60[ncen];
    TH1F *hEtaFake_20[2][ncen],*hEtaFake_30[2][ncen], *hEtaFake_40[2][ncen],*hEtaFake_45[2][ncen], *hEtaFake_50[2][ncen], *hEtaFake_60[2][ncen];
    TH1F *hPhiFake_20[2][ncen],*hPhiFake_30[2][ncen], *hPhiFake_40[2][ncen],*hPhiFake_45[2][ncen], *hPhiFake_50[2][ncen], *hPhiFake_60[2][ncen];  
    TH1F *hEtaFakeRatio_20[2][ncen],*hEtaFakeRatio_30[2][ncen], *hEtaFakeRatio_40[2][ncen],*hEtaFakeRatio_45[2][ncen], *hEtaFakeRatio_50[2][ncen], *hEtaFakeRatio_60[2][ncen];
    TH1F *hPhiFakeRatio_20[2][ncen],*hPhiFakeRatio_30[2][ncen], *hPhiFakeRatio_40[2][ncen],*hPhiFakeRatio_45[2][ncen], *hPhiFakeRatio_50[2][ncen], *hPhiFakeRatio_60[2][ncen];

    TH1F *hEtaFakeRatio_70[2][ncen], *hEtaFakeRatio_80[2][ncen];
    TH1F *hPhiFakeRatio_70[2][ncen], *hPhiFakeRatio_80[2][ncen];
    TH1F *hEtaFakeAll_70[ncen], *hEtaFakeAll_80[ncen];
    TH1F *hPhiFakeAll_70[ncen], *hPhiFakeAll_80[ncen];
    TH1F *hEtaFake_70[2][ncen], *hEtaFake_80[2][ncen];
    TH1F *hPhiFake_70[2][ncen], *hPhiFake_80[2][ncen];
    
  TFile *fin = new TFile(Form("Histos/OutputHist_ntuples_reduced_eta_lt1pt8_%s.root",kSpecies.c_str()),"r");
// radius loop
    for(int nj=0; nj<knj; nj++){
    // i is the with or without jetID loop, NOT the nj loop! 0 == no JetID, 1== JetID
     for(int i=0; i<2; i++){
 // centrality loop    
      for(int ic=0; ic<ncen; ic++){
    	if(i==0){
     	  hPtFakeAll[ic] = (TH1F*)fin->Get(Form("hPtFakeAll_%s_%d",(kAlgName+srad[nj]+kjetType).c_str(),ic)); 
     	  hPtFakeAll[ic]->Print("base");  	  
    	  hEtaFakeAll_20[ic] = (TH1F*)fin->Get(Form("hEtaFakeAll_20_%s_%d",(kAlgName+srad[nj]+kjetType).c_str(),ic)); 
    	  hEtaFakeAll_20[ic]->Print("base");
    	  hPhiFakeAll_20[ic] = (TH1F*)fin->Get(Form("hPhiFakeAll_20_%s_%d",(kAlgName+srad[nj]+kjetType).c_str(),ic)); 
          hPhiFakeAll_20[ic]->Print("base");
    	  hEtaFakeAll_30[ic] = (TH1F*)fin->Get(Form("hEtaFakeAll_30_%s_%d",(kAlgName+srad[nj]+kjetType).c_str(),ic)); 
    	  hEtaFakeAll_30[ic]->Print("base");
    	  hPhiFakeAll_30[ic] = (TH1F*)fin->Get(Form("hPhiFakeAll_30_%s_%d",(kAlgName+srad[nj]+kjetType).c_str(),ic)); 
          hPhiFakeAll_30[ic]->Print("base");
    	  hEtaFakeAll_40[ic] = (TH1F*)fin->Get(Form("hEtaFakeAll_40_%s_%d",(kAlgName+srad[nj]+kjetType).c_str(),ic)); 
    	  hEtaFakeAll_40[ic]->Print("base");
    	  hPhiFakeAll_40[ic] = (TH1F*)fin->Get(Form("hPhiFakeAll_40_%s_%d",(kAlgName+srad[nj]+kjetType).c_str(),ic)); 
          hPhiFakeAll_40[ic]->Print("base");
    	  hEtaFakeAll_45[ic] = (TH1F*)fin->Get(Form("hEtaFakeAll_45_%s_%d",(kAlgName+srad[nj]+kjetType).c_str(),ic)); 
    	  hEtaFakeAll_45[ic] ->Print("base");
    	  hPhiFakeAll_45[ic] = (TH1F*)fin->Get(Form("hPhiFakeAll_45_%s_%d",(kAlgName+srad[nj]+kjetType).c_str(),ic)); 
          hPhiFakeAll_45[ic]->Print("base");
    	  hEtaFakeAll_50[ic] = (TH1F*)fin->Get(Form("hEtaFakeAll_50_%s_%d",(kAlgName+srad[nj]+kjetType).c_str(),ic));
    	  hEtaFakeAll_50[ic]->Print("base"); 
    	  hPhiFakeAll_50[ic] = (TH1F*)fin->Get(Form("hPhiFakeAll_50_%s_%d",(kAlgName+srad[nj]+kjetType).c_str(),ic)); 
    	  hPhiFakeAll_50[ic]->Print("base"); 
    	  hEtaFakeAll_60[ic] = (TH1F*)fin->Get(Form("hEtaFakeAll_60_%s_%d",(kAlgName+srad[nj]+kjetType).c_str(),ic)); 
    	  hEtaFakeAll_60[ic]->Print("base");
    	  hPhiFakeAll_60[ic] = (TH1F*)fin->Get(Form("hPhiFakeAll_60_%s_%d",(kAlgName+srad[nj]+kjetType).c_str(),ic)); 
          hPhiFakeAll_60[ic]->Print("base");
    	  hEtaFakeAll_70[ic] = (TH1F*)fin->Get(Form("hEtaFakeAll_70_%s_%d",(kAlgName+srad[nj]+kjetType).c_str(),ic)); 
    	  hEtaFakeAll_70[ic]->Print("base");
    	  hPhiFakeAll_70[ic] = (TH1F*)fin->Get(Form("hPhiFakeAll_70_%s_%d",(kAlgName+srad[nj]+kjetType).c_str(),ic)); 
          hPhiFakeAll_70[ic]->Print("base");
    	  hEtaFakeAll_80[ic] = (TH1F*)fin->Get(Form("hEtaFakeAll_80_%s_%d",(kAlgName+srad[nj]+kjetType).c_str(),ic)); 
    	  hEtaFakeAll_80[ic]->Print("base");
    	  hPhiFakeAll_80[ic] = (TH1F*)fin->Get(Form("hPhiFakeAll_80_%s_%d",(kAlgName+srad[nj]+kjetType).c_str(),ic)); 
          hPhiFakeAll_80[ic]->Print("base");
           for(int in=0; in<neta; in++){    
              hPtFakeAll_etabin[ic][in] = (TH1F*)fin->Get(Form("hPtFakeAll_etabin_%s_%d_%d",(kAlgName+srad[nj]+kjetType).c_str(),ic,in)); 
     	      hPtFakeAll_etabin[ic][in]->Print("base");  
      	  }
   	  
    	}
     	  hPtAll[i][ic] = (TH1F*)fin->Get(Form("hPtAll_%s_%d_%d",(kAlgName+srad[nj]+kjetType).c_str(),i,ic));
     	  hPtAll[i][ic]->Print("base");    	
    	// end if i=0 (no jetID)
     	hPtEff[i][ic] = (TH1F*)fin->Get(Form("hPtEff_%s_%d_%d",(kAlgName+srad[nj]+kjetType).c_str(),i,ic)); 
     	hPtEff[i][ic]->Print("base");
    	hPtFake[i][ic] = (TH1F*)fin->Get(Form("hPtFake_%s_%d_%d",(kAlgName+srad[nj]+kjetType).c_str(),i,ic)); 
        hPtFake[i][ic]->Print("base");
    	hEtaFake_20[i][ic] = (TH1F*)fin->Get(Form("hEtaFake_20_%s_%d_%d",(kAlgName+srad[nj]+kjetType).c_str(),i,ic)); 
    	hEtaFake_20[i][ic]->Print("base");
    	hPhiFake_20[i][ic] = (TH1F*)fin->Get(Form("hPhiFake_20_%s_%d_%d",(kAlgName+srad[nj]+kjetType).c_str(),i,ic)); 
//.........这里部分代码省略.........

void global_reco(Int_t nEvents = 100, Int_t seed = 555)
{
	gRandom->SetSeed(seed);

   TTree::SetMaxTreeSize(90000000000);
   TString script = TString(gSystem->Getenv("LIT_SCRIPT"));
   TString parDir = TString(gSystem->Getenv("VMCWORKDIR")) + TString("/parameters");

   // Input and output data
   TString dir = "events/much_anna_omega_8gev_10k/"; // Output directory
   TString mcFile = dir + "mc.0000.root"; // MC transport file
   TString parFile = dir + "param.0000.root"; // Parameters file
   TString globalRecoFile = dir + "global.reco.0000.root"; // Output file with reconstructed tracks and hits

   // Digi files
   TList* parFileList = new TList();
   TObjString stsDigiFile = parDir + "/sts/sts_v12b_std.digi.par"; // STS digi file
   TObjString trdDigiFile = parDir + "/trd/trd_v13g.digi.par"; // TRD digi file
   TString muchDigiFile = parDir + "/much/much_v12c.digi.root"; // MUCH digi file
   TString stsMatBudgetFile = parDir + "/sts/sts_matbudget_v12b.root";
   TObjString tofDigiFile = parDir + "/tof/tof_v13b.digi.par";// TOF digi file

   // Reconstruction parameters
   TString globalTrackingType = "nn"; // Global tracking type
   TString trdHitProducerType = "smearing"; // TRD hit producer type: smearing, digi, clustering
   TString muchHitProducerType = "advanced"; // MUCH hit producer type: simple, advanced

   if (script == "yes") {
      mcFile = TString(gSystem->Getenv("LIT_MC_FILE"));
      parFile = TString(gSystem->Getenv("LIT_PAR_FILE"));
      globalRecoFile = TString(gSystem->Getenv("LIT_GLOBAL_RECO_FILE"));

      stsDigiFile = TString(gSystem->Getenv("LIT_STS_DIGI"));
      trdDigiFile = TString(gSystem->Getenv("LIT_TRD_DIGI"));
      muchDigiFile = TString(gSystem->Getenv("LIT_MUCH_DIGI"));
      tofDigiFile = TString(gSystem->Getenv("LIT_TOF_DIGI"));
      stsMatBudgetFile = TString(gSystem->Getenv("LIT_STS_MAT_BUDGET_FILE"));
   }

   parFileList->Add(&stsDigiFile);
   parFileList->Add(&trdDigiFile);
   parFileList->Add(&tofDigiFile);

	Int_t iVerbose = 1;
	TStopwatch timer;
	timer.Start();

	gROOT->LoadMacro("$VMCWORKDIR/macro/littrack/loadlibs.C");
	loadlibs();

	FairRunAna *run = new FairRunAna();
   run->SetInputFile(mcFile);
	run->SetOutputFile(globalRecoFile);


   // -----   STS  -------------------------------------------------
   Double_t threshold = 4;
   Double_t noiseWidth = 0.01;
   Int_t nofBits = 12;
   Double_t ElectronsPerAdc =  10.;
   Double_t StripDeadTime = 0.1;
   CbmStsDigitize* stsDigitize = new CbmStsDigitize("STS Digitiser", iVerbose);
   stsDigitize->SetRealisticResponse();
   stsDigitize->SetFrontThreshold(threshold);
   stsDigitize->SetBackThreshold(threshold);
   stsDigitize->SetFrontNoiseWidth(noiseWidth);
   stsDigitize->SetBackNoiseWidth(noiseWidth);
   stsDigitize->SetFrontNofBits(nofBits);
   stsDigitize->SetBackNofBits(nofBits);
   stsDigitize->SetFrontNofElPerAdc(ElectronsPerAdc);
   stsDigitize->SetBackNofElPerAdc(ElectronsPerAdc);
   stsDigitize->SetStripDeadTime(StripDeadTime);
   run->AddTask(stsDigitize);

   FairTask* stsClusterFinder = new CbmStsClusterFinder("STS Cluster Finder",iVerbose);
   run->AddTask(stsClusterFinder);

   FairTask* stsFindHits = new CbmStsFindHits("STS Hit Finder", iVerbose);
   run->AddTask(stsFindHits);

   FairTask* stsMatchHits = new CbmStsMatchHits("STS Hit Matcher", iVerbose);
   run->AddTask(stsMatchHits);

   FairTask* kalman = new CbmKF();
   run->AddTask(kalman);
   CbmL1* l1 = new CbmL1();
   //l1->SetExtrapolateToTheEndOfSTS(true);
   l1->SetMaterialBudgetFileName(stsMatBudgetFile);
   run->AddTask(l1);
   CbmStsTrackFinder* trackFinder = new CbmL1StsTrackFinder();
   FairTask* findTracks = new CbmStsFindTracks(iVerbose, trackFinder);
   run->AddTask(findTracks);

   FairTask* stsMatchTracks = new CbmStsMatchTracks("STSMatchTracks", iVerbose);
   run->AddTask(stsMatchTracks);
	// ------------------------------------------------------------------------

   if (IsMuch(parFile)) {
      // -------- MUCH digitization ------------
      CbmMuchDigitizeGem* digitize = new CbmMuchDigitizeGem(muchDigiFile.Data());
//.........这里部分代码省略.........

void chi2_lambda_stack_mumue(){
  TStopwatch totaltime;
  totaltime.Start();
  int sample, isPU, ischi2;
  
  TChain *chain_mu[N_chain];
  TChain *chain_el[N_chain];
  TChain *chain_jet[N_chain];
  TChain *chain_event[N_chain];
  
  cout << "==========================" << endl
  << "40, 50, 60, 100(no gen lvl info) : signal" << endl
  << "==========================" << endl
  << "Run : ";
  cin >> sample;
  cout << endl;
  
  cout << "==========================" << endl
  << "PileUP? (0 = no, 1 = yes) : ";
  cin >> isPU;
  cout << "==========================" << endl
  << "chi2 fit? (0 = no, 1 = yes) :";
  cin >> ischi2;
  cout << endl;

  /*
  /// [0] = Signal ///
  chain_mu[0] = new TChain("Muon");
  chain_el[0] = new TChain("Electron");
  chain_event[0] = new TChain("Event");
  chain_mu[0]->Add("./files/signal/trilepton_mumue_SKHN"+TString::Itoa(sample, 10)+"_lep_5_3_14.root");
  chain_el[0]->Add("./files/signal/trilepton_mumue_SKHN"+TString::Itoa(sample, 10)+"_lep_5_3_14.root");
  chain_event[0]->Add("./files/signal/trilepton_mumue_SKHN"+TString::Itoa(sample, 10)+"_lep_5_3_14.root");
  */
  
  /// [1] = Data ///
  chain_mu[1] = new TChain("Muon");
  chain_el[1] = new TChain("Electron");
  chain_event[1] = new TChain("Event");
  chain_jet[1] = new TChain("Jet");
  chain_mu[1]->Add("./files/data/*.root");
  chain_el[1]->Add("./files/data/*.root");
  chain_event[1]->Add("./files/data/*.root");
  chain_jet[1]->Add("./files/data/*.root");
  
  /// [2]~ = Bkg ///
  TString chainlist[N_chain] = {"signal", "data",
                                "DY10to50", "DY50plus", "Zbb",  // DY
                                "WZtollln_mg", "WZtollqq_mg", "WZtoqqln_mg", "ZZtollll_mg", "ZZtollnn_mg", "ZZtollqq_mg", "WW_mg",  // VV
                                "Wbb", // Wjets
                                "topDIL", // Top
                                "TTG", "TTWW", "WWG", "WWW", "WWZ", "WZZ", "ZZZ", "ttZ", "HtoWW", "Wtollln"}; // others
  for(int i=2; i<N_chain; i++){
    chain_mu[i] = new TChain("Muon");
    chain_el[i] = new TChain("Electron");
    chain_event[i] = new TChain("Event");
    chain_jet[i] = new TChain("Jet");
    chain_mu[i]->Add("./files/bkg/trilepton_mumue_SK"+chainlist[i]+"_5_3_14.root");
    chain_el[i]->Add("./files/bkg/trilepton_mumue_SK"+chainlist[i]+"_5_3_14.root");
    chain_event[i]->Add("./files/bkg/trilepton_mumue_SK"+chainlist[i]+"_5_3_14.root");
    chain_jet[i]->Add("./files/bkg/trilepton_mumue_SK"+chainlist[i]+"_5_3_14.root");
  }
  
  Muon *mu[N_chain];
  Electron *el[N_chain];
  Event *evt[N_chain];
  Jet *jet[N_chain];
           
  double MET, METphi, nu_Pz, W_mass_fit, MET_fit, weight, PUweight_sum;
  int N_entry, HN_x_min=0, HN_x_max=100, HN_dx=10, W_x_min=70, W_x_max=120, W_dx=5, solution_selection, smaller;
  struct Fit fitresult[2];
  Color_t fillcolors[N_chain] = {0, kRed-5, // signal and data. data random color
    kAzure-5, kAzure-5, kAzure-5,
    kGreen-5, kGreen-5, kGreen-5, kGreen-5, kGreen-5, kGreen-5, kGreen-5,
    kYellow-5,
    kRed-5,
    kGreen-10, kGreen-10, kGreen-10, kGreen-10, kGreen-10, kGreen-10, kGreen-10, kGreen-10, kGreen-10, kOrange-4};
  double gamma_star_bin[16];
  TString gamma_star_bin_label[15];
  for(int i=0; i<15; i++){
    gamma_star_bin[i] = 0.5*i;
    if(i%2==0) gamma_star_bin_label[i] = "";
    else      gamma_star_bin_label[i] = "            "+TString::Itoa(i/2+1, 10)+".0";
    //cout << gamma_star_bin[i] << '\t' << gamma_star_bin_label[i] << endl;
  }
  gamma_star_bin[15] = 8;
  gamma_star_bin_label[14] = "        > 7.0 GeV";
  
  TH1D *hist_dimuon[N_chain];
  
  TH1D *hist_dimuon_err = new TH1D("hist_dimuon_err", "", (50.-0.)/1., 0., 50.);
  
  THStack *bkgstack_dimuon = new THStack("bkgstack_dimuon", "");
  
  TLorentzVector mu_4vec[2], el_4vec, W_real, HN, nu, gen_nu, gamma_star, z_candidate, selection_nu[2];

  for(int k=1; k<N_chain; k++){
    N_entry = 0;
    PUweight_sum=0;
    
//.........这里部分代码省略.........

//.........这里部分代码省略.........
std::vector<double> *myweight=0;
std::vector<double> *myjetmatchinginfo=0;
std::vector<double> *mycoswphoton=0;

   std::cout << "--- Select signal sample" << std::endl;
   TTree* theTree = (TTree*)input->Get("analyzestep2/atq");
//   Int_t myjetmultiplicity, mybjetmultiplicity , myleptoncharge;
//   Float_t userVar1, userVar2;

   theTree->SetBranchAddress("ptphoton", &myptphoton  );
   theTree->SetBranchAddress( "etaphoton", &myetaphoton );
   theTree->SetBranchAddress( "ptmuon", &myptmuon );
   theTree->SetBranchAddress( "etamuon", &myetamuon );
   theTree->SetBranchAddress( "ptjet", &myptjet );
   theTree->SetBranchAddress( "etajet", &myetajet );
   theTree->SetBranchAddress( "masstop", &mymasstop );
//   theTree->SetBranchAddress( "mtw", &mymtw );
   theTree->SetBranchAddress( "deltaRphotonjet", &mydeltaRphotonjet );
   theTree->SetBranchAddress( "deltaRphotonmuon", &mydeltaRphotonmuon );
//   theTree->SetBranchAddress( "ht", &myht );
   theTree->SetBranchAddress( "costopphoton", &mycostopphoton );
   theTree->SetBranchAddress( "jetmultiplicity", &myjetmultiplicity );
//   theTree->SetBranchAddress( "bjetmultiplicity", &mybjetmultiplicity );
   theTree->SetBranchAddress( "deltaphiphotonmet", &mydeltaphiphotonmet );
   theTree->SetBranchAddress( "cvsdiscriminant", &mycvsdiscriminant );
//   theTree->SetBranchAddress( "leptoncharge", &myleptoncharge );
   theTree->SetBranchAddress( "weight", &myweight);
   theTree->SetBranchAddress( "coswphoton", &mycoswphoton );
	theTree->SetBranchAddress( "jetmatchinginfo", &myjetmatchinginfo );


//   std::cout << "--- Processing: " << theTree->GetEntries() << " events" << std::endl;
   TStopwatch sw;
   sw.Start();
   for (Long64_t ievt=0; ievt<theTree->GetEntries();ievt++) {
//   std::cout << "--- ... Processing event: " << ievt << std::endl;
double finalweight;

      if (ievt%1000 == 0) std::cout << "--- ... Processing event: " << ievt << std::endl;

      theTree->GetEntry(ievt);
//for (int l=0;l<sizeof(myptphoton);l++){
//std::cout << "--- ... reza: " << myptphoton[l] <<std::endl;
//}
//std::cout << "--- ......................."<< (*mycvsdiscriminant)[0]<<std::endl;
      // --- Return the MVA outputs and fill into histograms

finalweight=(*myweight)[0];
//cout<<(*myweight)[0]<<endl;
if((*mymasstop )[0]>mtopdown && (*mymasstop )[0]<mtopup){
hists[idx] ->Fill( (*myetajet)[0],finalweight );
if (samples_[idx]=="WPHJET.root")insidewphjet=insidewphjet+finalweight;
if (samples_[idx]=="SIGNALtGu.root")nsignalevent=nsignalevent+1;
//cout<<insidewphjet<<endl;
}
else {
SBhists[idx] ->Fill( (*myetajet)[0],finalweight );
if (samples_[idx]=="WPHJET.root")outsidewphjet=outsidewphjet+finalweight;}


      // Retrieve also per-event error
}
delete myptphoton;
delete myetaphoton;
delete myptmuon;
delete myetamuon;

//.........这里部分代码省略.........
  
   // run first a data fit 
  
   const RooArgSet * poiSet = sbModel->GetParametersOfInterest();
   RooRealVar *poi = (RooRealVar*)poiSet->first();
  
   std::cout << "StandardHypoTestInvDemo : POI initial value:   " << poi->GetName() << " = " << poi->getVal()   << std::endl;  
  
   // fit the data first (need to use constraint )
   TStopwatch tw; 

   bool doFit = initialFit;
   if (testStatType == 0 && initialFit == -1) doFit = false;  // case of LEP test statistic
   if (type == 3  && initialFit == -1) doFit = false;         // case of Asymptoticcalculator with nominal Asimov
   double poihat = 0;

   if (minimizerType.size()==0) minimizerType = ROOT::Math::MinimizerOptions::DefaultMinimizerType();
   else 
      ROOT::Math::MinimizerOptions::SetDefaultMinimizer(minimizerType.c_str());
    
   Info("StandardHypoTestInvDemo","Using %s as minimizer for computing the test statistic",
        ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str() );
   
   if (doFit)  { 

      // do the fit : By doing a fit the POI snapshot (for S+B)  is set to the fit value
      // and the nuisance parameters nominal values will be set to the fit value. 
      // This is relevant when using LEP test statistics

      Info( "StandardHypoTestInvDemo"," Doing a first fit to the observed data ");
      RooArgSet constrainParams;
      if (sbModel->GetNuisanceParameters() ) constrainParams.add(*sbModel->GetNuisanceParameters());
      RooStats::RemoveConstantParameters(&constrainParams);
      tw.Start(); 
      RooFitResult * fitres = sbModel->GetPdf()->fitTo(*data,InitialHesse(false), Hesse(false),
                                                       Minimizer(minimizerType.c_str(),"Migrad"), Strategy(0), PrintLevel(mPrintLevel), Constrain(constrainParams), Save(true) );
      if (fitres->status() != 0) { 
         Warning("StandardHypoTestInvDemo","Fit to the model failed - try with strategy 1 and perform first an Hesse computation");
         fitres = sbModel->GetPdf()->fitTo(*data,InitialHesse(true), Hesse(false),Minimizer(minimizerType.c_str(),"Migrad"), Strategy(1), PrintLevel(mPrintLevel+1), Constrain(constrainParams), Save(true) );
      }
      if (fitres->status() != 0) 
         Warning("StandardHypoTestInvDemo"," Fit still failed - continue anyway.....");
  
  
      poihat  = poi->getVal();
      std::cout << "StandardHypoTestInvDemo - Best Fit value : " << poi->GetName() << " = "  
                << poihat << " +/- " << poi->getError() << std::endl;
      std::cout << "Time for fitting : "; tw.Print(); 
  
      //save best fit value in the poi snapshot 
      sbModel->SetSnapshot(*sbModel->GetParametersOfInterest());
      std::cout << "StandardHypoTestInvo: snapshot of S+B Model " << sbModel->GetName() 
                << " is set to the best fit value" << std::endl;
  
   }

   // print a message in case of LEP test statistics because it affects result by doing or not doing a fit 
   if (testStatType == 0) {
      if (!doFit) 
         Info("StandardHypoTestInvDemo","Using LEP test statistic - an initial fit is not done and the TS will use the nuisances at the model value");
      else 
         Info("StandardHypoTestInvDemo","Using LEP test statistic - an initial fit has been done and the TS will use the nuisances at the best fit value");
   }


   // build test statistics and hypotest calculators for running the inverter 

void ana_standalone_single_photon_MC() {

  gSystem->Load("/afs/cern.ch/user/y/yohay/scratch0/CMSSW_4_2_4_patch2/src/SusyAnalysis/SusyNtuplizer/macro/libSusy.so");

  // // Look ../jec/JetMETObjects/README
  gSystem->Load("/afs/cern.ch/user/y/yohay/scratch0/CMSSW_4_2_4_patch2/src/SusyAnalysis/SusyNtuplizer/jec/lib/libJetMETObjects.so");

  //GMSBTools shared library (contains Categorizer class)
  gSystem->Load("/afs/cern.ch/user/y/yohay/scratch0/CMSSW_4_2_4_patch2/src/GMSBTools/lib/libFilters.so");

  // Printing utility for ntuple variables
  gSystem->Load("/afs/cern.ch/user/y/yohay/scratch0/CMSSW_4_2_4_patch2/src/SusyAnalysis/SusyNtuplizer/macro/SusyEventPrinter_cc.so");

  // Main analysis code
  gSystem->SetIncludePath("-I/afs/cern.ch/user/y/yohay/scratch0/CMSSW_4_2_4_patch2/src");
  gSystem->Load("/afs/cern.ch/user/y/yohay/scratch0/CMSSW_4_2_4_patch2/src/SusyAnalysis/SusyNtuplizer/macro/EventAnalyzer_cc.so");
  gSystem->Load("/afs/cern.ch/user/y/yohay/scratch0/CMSSW_4_2_4_patch2/src/SusyAnalysis/SusyNtuplizer/macro/DongwookCategoryProducer_cc.so");

  /*now, ECAL/HCAL/track isolation cuts are meant to be matched with loose trigger cuts, so don't 
    use them (i.e. track isolation) to distinguish photon from fake!*/

  //configuration
  ParameterSet pars;
  pars.photonTag = "photons";
  pars.photon1ETMin = 40.0;
  pars.photon2ETMin = 30.0;
  pars.photonAbsEtaMax = 1.4442;
  pars.photonECALIsoMaxPTMultiplier = 0.012; /*https://twiki.cern.ch/twiki/bin/viewauth/CMS/
					       EgammaWorkingPointsv3, IsoVL, assuming dR = 0.3 
					       cone*/
  pars.photonECALIsoMaxConstant = 6.0; /*https://twiki.cern.ch/twiki/bin/viewauth/CMS/
					 EgammaWorkingPointsv3, IsoVL, assuming dR = 0.3 cone*/
  pars.photonECALIsoEffArea = 0.1474; /*https://twiki.cern.ch/twiki/bin/view/CMS/
					RA3IsolationConePileupCorrections, dR = 0.3 cone, 
					18-Oct-11*/
  pars.photonHCALIsoMaxPTMultiplier = 0.005; /*https://twiki.cern.ch/twiki/bin/viewauth/CMS/
					       EgammaWorkingPointsv3, IsoVL, assuming dR = 0.3 
					       cone*/
  pars.photonHCALIsoMaxConstant = 4.0; /*https://twiki.cern.ch/twiki/bin/viewauth/CMS/
					 EgammaWorkingPointsv3, IsoVL, assuming dR = 0.3 cone*/
  pars.photonHCALIsoEffArea = 0.0467; /*https://twiki.cern.ch/twiki/bin/view/CMS/
					RA3IsolationConePileupCorrections, dR = 0.3 cone, 
					18-Oct-11*/
  pars.photonHOverEMax = 0.05;
  pars.photonR9Max = 0.98;
  pars.photonTrackIsoMaxPTMultiplier = 0.002; /*https://twiki.cern.ch/twiki/bin/viewauth/CMS/
						EgammaWorkingPointsv3, IsoVL, assuming dR = 0.3 
						cone*/
  pars.photonTrackIsoMaxConstant = 4.0; /*https://twiki.cern.ch/twiki/bin/viewauth/CMS/
					  EgammaWorkingPointsv3, IsoVL, assuming dR = 0.3 cone*/
  pars.photonCombinedIsoMax = 6.0;
  pars.fakeCombinedIsoMax = 12.0;
  pars.isoConeHLT = DR03;
  pars.isoConeOffline = DR03;
  pars.photonSigmaIetaIetaMax = 0.011;
  pars.photonAbsSeedTimeMax = -1.0;
  pars.photonE2OverE9Max = -1.0;
  pars.photonDPhiMin = 0.05;
  pars.photonDRMin = 0.8;
  pars.pixelVetoOnFake = true;
  pars.treeName = "susyTree";
//   pars.input = VSTRING(1, "/data2/yohay/RA3/Data2011A_ToRun167913_Filter-JsonHLTtwo43-30GeVPhosWithR9HoverE_NoPileupCorr_Photon_NEW.root");
  pars.input = VSTRING();
FILES
  pars.HLT = vector<TString>();
  pars.HLT.push_back("HLT_Photon70_CaloIdL_HT300");
  pars.HLT.push_back("HLT_Photon70_CaloIdL_HT400");
  pars.HLT.push_back("HLT_Photon135");
  pars.HLT.push_back("HLT_Photon75_CaloIdVL_IsoL");
  pars.HLT.push_back("HLT_Photon125");
  pars.nEvts = -1;
  pars.JSON = "";
//   pars.JSON = "/afs/cern.ch/user/y/yohay/scratch0/CMSSW_4_2_4_patch2/src/SusyAnalysis/SusyNtuplizer/macro/JSON_160431-177878_May10ReReco_Run2011APromptRecov4v6_Aug5ReReco_Run2011BPromptRecov1.txt";
//   pars.outputFile = "/data2/yohay/RA3/Data2011A_ToRun167913_Filter-JsonHLTtwo43-30GeVPhosWithR9HoverE_NoPileupCorr_Photon_NEW_categorized_OR.root";
//   pars.outputFile = "/data2/yohay/RA3/1140pb-1_ff_categorized_new.root";
  pars.outputFile = "DATASET_JSON_HLT_PV_single_photon_skim.root";

  TStopwatch ts;

  ts.Start();

  DongwookCategoryProducer producer(pars);

  ts.Stop();

  std::cout << "RealTime : " << ts.RealTime()/60.0 << " minutes" << std::endl;
  std::cout << "CPUTime  : " << ts.CpuTime()/60.0 << " minutes" << std::endl;

}

void 
Frequentist_calculation(const char * infile = "results/tttt_combined_my_measurement_model.root",
                        const char * wsName = "combined",
                        const char * modelSBName = "ModelConfig",
                        const char * modelBName = "",
                        const char * dataName = "obsData",                 
                        int calculatorType = 0,
                        int testStatType = 3, 
                        bool useCLs = true ,  
                        int npoints = 25,
                        double poimin = 0.,
                        double poimax = 100.,
                        int ntoys=1000,
                        bool useNumberCounting = false,
                        const char * nuisPriorName = 0){
/*

  Other Parameter to pass in tutorial
  apart from standard for filename, ws, modelconfig and data

  (calculator)Type = 0 Freq calculator 
  (calculator)Type = 1 Hybrid calculator
  (calculator)Type = 2 Asymptotic calculator  
  (calculator)Type = 3 Asymptotic calculator using nominal Asimov data sets (not using fitted parameter values but nominal ones)

  testStatType = 0 LEP
  = 1 Tevatron 
  = 2 Profile Likelihood
  = 3 Profile Likelihood one sided (i.e. = 0 if mu < mu_hat)
  = 4 Profiel Likelihood signed ( pll = -pll if mu < mu_hat) 
  = 5 Max Likelihood Estimate as test statistic
  = 6 Number of observed event as test statistic

  useCLs          scan for CLs (otherwise for CLs+b)    

  npoints:        number of points to scan , for autoscan set npoints = -1 

  poimin,poimax:  min/max value to scan in case of fixed scans 
  (if min >  max, try to find automatically)                           

  ntoys:         number of toys to use 

  useNumberCounting:  set to true when using number counting events 

  nuisPriorName:   name of prior for the nnuisance. This is often expressed as constraint term in the global model
  It is needed only when using the HybridCalculator (type=1)
  If not given by default the prior pdf from ModelConfig is used. 

  extra options are available as global paramwters of the macro. They major ones are: 

  plotHypoTestResult   plot result of tests at each point (TS distributions) (defauly is true)
  useProof             use Proof   (default is true) 
  writeResult          write result of scan (default is true)
  rebuild              rebuild scan for expected limits (require extra toys) (default is false)
  generateBinned       generate binned data sets for toys (default is false) - be careful not to activate with 
  a too large (>=3) number of observables 
  nToyRatio            ratio of S+B/B toys (default is 2)
    

*/


  // Time this macro
  TStopwatch t;
  t.Start();
  
  
   TString fileName(infile);
   if (fileName.IsNull()) { 
      fileName = "results/example_combined_GaussExample_model.root";
      std::cout << "Use standard file generated with HistFactory : " << fileName << std::endl;
   }
  
   // open file and check if input file exists
   TFile * file = TFile::Open(fileName); 
  
   // if input file was specified but not found, quit
   if(!file && !TString(infile).IsNull()){
      cout <<"file " << fileName << " not found" << endl;
      return;
   } 
  
/*   // if default file not found, try to create it
   if(!file ){
      // Normally this would be run on the command line
      cout <<"will run standard hist2workspace example"<<endl;
      gROOT->ProcessLine(".! prepareHistFactory .");
      gROOT->ProcessLine(".! hist2workspace config/example.xml");
      cout <<"\n\n---------------------"<<endl;
      cout <<"Done creating example input"<<endl;
      cout <<"---------------------\n\n"<<endl;
    
      // now try to access the file again
      file = TFile::Open(fileName);
    
   }
*/  
   if(!file){
      // if it is still not there, then we can't continue
      cout << "Not able to run hist2workspace to create example input" <<endl;
//.........这里部分代码省略.........

void r3ball_batch(Int_t nEvents = 1,
	    TObjArray& fDetList,
	    TString Target = "LeadTarget",
		Bool_t fVis=kFALSE,
		TString fMC="TGeant3",
	    TString fGenerator="box",
	    Bool_t fUserPList= kFALSE,
		Bool_t fR3BMagnet= kTRUE,
		Double_t fEnergyP=1.0,
		Int_t fMult=1, 
		Int_t fGeoVer=5, 
		Double_t fNonUni=1.0	
			)
{


  TString dir = getenv("VMCWORKDIR");
  TString r3bdir = dir + "/macros";

  TString r3b_geomdir = dir + "/geometry";
  gSystem->Setenv("GEOMPATH",r3b_geomdir.Data());

  TString r3b_confdir = dir + "gconfig";
  gSystem->Setenv("CONFIG_DIR",r3b_confdir.Data());

// Output files
  TString OutFile = "r3bsim.root";
  TString ParFile = "r3bpar.root";


  // In general, the following parts need not be touched
  // ========================================================================

  // ----    Debug option   -------------------------------------------------
  gDebug = 0;
  // ------------------------------------------------------------------------

  // -----   Timer   --------------------------------------------------------
  TStopwatch timer;
  timer.Start();
  // ------------------------------------------------------------------------

 
  // -----   Create simulation run   ----------------------------------------
  FairRunSim* run = new FairRunSim();
  run->SetName(fMC.Data());                  // Transport engine
  run->SetOutputFile(OutFile.Data());        // Output file
  FairRuntimeDb* rtdb = run->GetRuntimeDb();

  //  R3B Special Physics List in G4 case
  if ( (fUserPList  == kTRUE ) &&
       (fMC.CompareTo("TGeant4")   == 0)
      ){
       run->SetUserConfig("g4R3bConfig.C");
       run->SetUserCuts("SetR3BCuts.C");
   }


  // -----   Create media   -------------------------------------------------
  run->SetMaterials("media_r3b.geo");       // Materials
  

  // Magnetic field map type
   Int_t fFieldMap = 0;


  // Global Transformations
  //- Two ways for a Volume Rotation are supported
  //-- 1) Global Rotation (Euler Angles definition)
  //-- This represent the composition of : first a rotation about Z axis with
  //-- angle phi, then a rotation with theta about the rotated X axis, and
  //-- finally a rotation with psi about the new Z axis.
  Double_t phi,theta,psi;

  //-- 2) Rotation in Ref. Frame of the Volume
  //-- Rotation is Using Local Ref. Frame axis angles
  Double_t thetaX,thetaY,thetaZ;


  //- Global Translation  Lab. frame.
  Double_t tx,ty,tz;

  // - Polar angular limits
  Double_t minTheta=35., maxTheta=55.;


  // -----   Create R3B geometry --------------------------------------------
  //R3B Cave definition
  FairModule* cave= new R3BCave("CAVE");
  cave->SetGeometryFileName("r3b_cave.geo");
  run->AddModule(cave);


  //R3B Target definition
  if (fDetList.FindObject("TARGET") ) {
      R3BModule* target= new R3BTarget(Target.Data());
      // Global Lab. Rotation
      phi    =  0.0; // (deg)
      theta  =  0.0; // (deg)
      psi    =  0.0; // (deg)
//.........这里部分代码省略.........

//______________________________________________________________________________
void AODtrainsim(Int_t merge=0)
{
  // Main analysis train macro.
  // merge = 0: production
  // merge = 1: intermediate merging
  // merge = 2: final merging + terminate
  // merge = 3: terminate only

  if (merge) {
    TGrid::Connect("alien://");
    if (!gGrid || !gGrid->IsConnected()) {
      ::Error("AODtrainsim.C", "No grid connection");
      return;
    }
  }
  // Set temporary merging directory to current one
  gSystem->Setenv("TMPDIR", gSystem->pwd());
  // Set temporary compilation directory to current one
  gSystem->SetBuildDir(gSystem->pwd(), kTRUE);
  printf("==================================================================\n");
  printf("===========    RUNNING FILTERING TRAIN   ==========\n");
  printf("==================================================================\n");
  printf("=  Configuring analysis train for:                               =\n");
  if (usePhysicsSelection)   printf("=  Physics selection                                                =\n");
  if (iESDfilter)   printf("=  ESD filter                                                    =\n");
  if (iMUONcopyAOD) printf("=  MUON copy AOD                                                 =\n");

  // Load common libraries and set include path
  gSystem->AddIncludePath("-I$ALICE_ROOT/include -I$ALICE_PHYSICS/include");

  // Make the analysis manager and connect event handlers
  AliAnalysisManager *mgr  = new AliAnalysisManager("Analysis Train", "Production train");
  if (useSysInfo) mgr->SetNSysInfo(100);

  // Load ParFiles
  if (!LoadAnalysisLibraries()) {
    ::Error("AODtrainsim.C", "Could not load analysis libraries");
    return;
  }
  
  // Create input handler (input container created automatically)
  // ESD input handler
  AliESDInputHandler *esdHandler = new AliESDInputHandler();
  mgr->SetInputEventHandler(esdHandler);
  // Monte Carlo handler
  if (useMC) {
    AliMCEventHandler* mcHandler = new AliMCEventHandler();
    mgr->SetMCtruthEventHandler(mcHandler);
    mcHandler->SetReadTR(useTR);
  }
  // AOD output container, created automatically when setting an AOD handler
  if (iAODhandler) {
    // AOD output handler
    AliAODHandler* aodHandler   = new AliAODHandler();
    aodHandler->SetOutputFileName("AliAOD.root");
    mgr->SetOutputEventHandler(aodHandler);
  }
  // Debugging if needed
  if (useDBG) mgr->SetDebugLevel(3);

  AddAnalysisTasks(merge);
  if (merge) {
    if (merge < 3) AODmerge();
    if (merge > 1) {
      mgr->InitAnalysis();
      mgr->SetGridHandler(new AliAnalysisAlien());
      mgr->StartAnalysis("grid terminate",0);
    }
    return;
  }
  // Run the analysis
  //
  TChain *chain = CreateChain();
  if (!chain) return;

  TStopwatch timer;
  timer.Start();
  mgr->SetSkipTerminate(kTRUE);
  if (mgr->InitAnalysis()) {
    mgr->PrintStatus();
    mgr->StartAnalysis("local", chain);
  }
  timer.Print();
}

void lhsig0810()
{

  TStopwatch timer;
  timer.Start();
  // Define histogarams
  
  gStyle->SetOptFit(1111);
  gStyle->SetOptStat(111111);
  
  TGraphErrors *g1 = new TGraphErrors(20);
  g1->SetName("Dilvsdm");
  g1->SetTitle("Dilution vs. dm");
  
  TGraph *g2 = new TGraph(20);
  g2->SetName("DilErrvsdm");
  g2->SetTitle("Dil Err vs. dm");
  
  for(Int_t i=0.;i<20.;i++){


    oscpar1_init = 1+i;

    cout << "============== dm fixed at "<< oscpar1_init <<" ===================="<<endl;

    //TH1F *h1 = new TH1F("h1","Gaussian Dist",50, 4.8, 5.8);
    TH1F *h2 = new TH1F("h2","Lifetime Dist",nBins, min, max);
    TH1F *h3 = new TH1F("h3","Lifetime Dist tag=1",nBins, min, max);
    TH1F *h4 = new TH1F("h4","Lifetime Dist tag=-1",nBins, min, max);
    
    // Get data
    
    // Generate events
    
    mixmasta_mc();
    
    for (Int_t ja=0; ja<nEvts; ja++){
      h2->Fill(lifetime[ja]);
      if (tag[ja] == 1){
	h3->Fill(lifetime[ja]);
      }else if(tag[ja] == -1){
	h4->Fill(lifetime[ja]);
      }else{
	cout << "Tag value "<<tag[ja]<< " out of range, should be -1 or 1" << endl;
	break;
      }
    }
    
    // Do unbinned likelihood fit
    
    TF1 *f3 = new TF1("f3", lftmosc_plt_d, min, max, 5);
    TF1 *f4 = new TF1("f4", lftmosc_plt_d, min, max, 5);
    
    

    
    unbinFitosc_d();
    
    for (Int_t j=0; j<4; j++){
      f3->SetParameter(j,fitpar[j]);
      f4->SetParameter(j,fitpar[j]);
    }
    f3->SetParameter(4,1);
    f4->SetParameter(4,-1);
    
    g1->SetPoint(i,fitpar[1],-(1-2*fitpar[3])/dilfit);
    g1->SetPointError(i,fiterr[1],2*fiterr[3]/dilfit);
    
    g2->SetPoint(i,fitpar[1],1.65*2*fiterr[3]/dilfit);
    
    delete h2;
    delete h3;
    delete h4;
    delete f3;
    delete f4;
  }
  
  TCanvas *vardm = new TCanvas("vardm","lhsig0810",800,400);
  vardm->Divide(2,1);
  
  vardm->cd(1);
  g1->GetXaxis()->SetTitle("dm");
  g1->GetXaxis()->CenterTitle();
  g1->GetYaxis()->SetTitle("Dilution (1-2alpha)");
  g1->GetYaxis()->CenterTitle();
  //g1->SetMarkerStyle(21);
  //g1->SetMarkerSize(1);
  g1->Draw("AP*"); 
  
  
  vardm->cd(2);
  
  gStyle->SetPadColor(10);
  gStyle->SetCanvasColor(10);
  vardm->SetGrid();
  
  g2->GetXaxis()->SetTitle("dm");
  g2->GetXaxis()->CenterTitle();
  g2->GetYaxis()->SetTitle("Dil. Err");
  g2->GetYaxis()->CenterTitle();
//.........这里部分代码省略.........

void call_validation()
{
  TStopwatch* myWatch = new TStopwatch();
  myWatch->Start();

//   std::string inFile = "DoubleElectron_Run2012A_filteredtree.root";
  
//   std::string inFile = "GluGluToHToZZTo2L2Q_M-400_8TeV-powheg-pythia6_filteredtree.root";
  std::string inFile = "DoubleMu_Run2012A_filteredtree.root";

  /*
  displayTwoTreesRatio(inFile.data(),"genParM_",100,0,1000,false,"test",
  		       "genParId_==25");
   displayTwoTreesRatio(inFile.data(),"PU_nTrueInt",60,0,60);

  displayTwoTreesRatio(inFile.data(),"EvtInfo_VertexX_[0]",100,0.23,0.26,false,"EvtInfo_VertexX_");
  displayTwoTreesRatio(inFile.data(),"EvtInfo_VertexY_[0]",100,0.38,0.41,false,"EvtInfo_VertexY_");

//   displayTwoTreesRatio(inFile.data(),"EvtInfo_VertexX_[0]",100,0.05,0.10,false,
// 		       "EvtInfo_VertexX_");
//   displayTwoTreesRatio(inFile.data(),"EvtInfo_VertexY_[0]",100,0.04,0.08,false,
// 		       "EvtInfo_VertexY_");
  
  displayTwoTreesRatio(inFile.data(),"EvtInfo_EventNum",100,0,100000000);
  displayTwoTreesRatio(inFile.data(),"EvtInfo_RunNum",100,150000,300000);
  displayTwoTreesRatio(inFile.data(),"EvtInfo_VertexZ_[0]",100,-30,30,false,
  "EvtInfo_VertexZ_");
  */
  displayTwoTreesRatio(inFile.data(),"EvtInfo_NumVtx",30,0,30);

  /*
  displayTwoTreesRatio(inFile.data(),"patJetPfAk05Pt_",50,0, 500);
  displayTwoTreesRatio(inFile.data(),"patJetPfAk05Eta_",50,-3.0,3.0);
  displayTwoTreesRatio(inFile.data(),"patJetPfAk05Phi_",50,0,TMath::Pi());
  displayTwoTreesRatio(inFile.data(),"patJetPfAk05CharMulti_",100,0,100);
  displayTwoTreesRatio(inFile.data(),"eleRho",100,0,70);
  displayTwoTreesRatio(inFile.data(),"muoRho",100,0,10);

  displayTwoTreesRatio(inFile.data(),"higgsM",50,0,1500);
  displayTwoTreesRatio(inFile.data(),"higgsPt",50, 0,500);
  displayTwoTreesRatio(inFile.data(),"higgsEta",50, -6,6);
  displayTwoTreesRatio(inFile.data(),"higgsPhi",50, 0,TMath::Pi());

  displayTwoTreesRatio(inFile.data(),"zllPt",50, 0,500);
  displayTwoTreesRatio(inFile.data(),"zllEta",50, -6,6);
  displayTwoTreesRatio(inFile.data(),"zllPhi",50, 0,TMath::Pi());

  displayTwoTreesRatio(inFile.data(),"zjjPt",50, 0,500);
  displayTwoTreesRatio(inFile.data(),"zjjEta",50, -6,6);
  displayTwoTreesRatio(inFile.data(),"zjjPhi",50, 0,TMath::Pi());

  displayTwoTreesRatio(inFile.data(),"zllM",100,60,120);
  displayTwoTreesRatio(inFile.data(),"zjjM",100, 0,300);

  displayTwoTreesRatio(inFile.data(),"nBTags",3,-0.5,2.5);
  displayTwoTreesRatio(inFile.data(),"lepType",2,-0.5,1.5);
  displayTwoTreesRatio(inFile.data(),"heliLD",100,0.0,1.0);
  */
  displayTwoTreesRatio(inFile.data(),"metSig", 50,0, 50,true);
  /*
  displayTwoTreesRatio(inFile.data(),"patJetPfAk05NeutEmEFr_",100,0,1,true);
  displayTwoTreesRatio(inFile.data(),"patJetPfAk05NeutHadEFr_",100,0,1,true);
  displayTwoTreesRatio(inFile.data(),"patJetPfAk05CharEmEFr_",100,0,1,true);
  displayTwoTreesRatio(inFile.data(),"patJetPfAk05CharHadEFr_",100,0,1,true);

  displayTwoTreesRatio(inFile.data(),"jetPt",50,0, 500);
  displayTwoTreesRatio(inFile.data(),"jetEta",50,-3.0,3.0);
  displayTwoTreesRatio(inFile.data(),"jetPhi",50,0,TMath::Pi());
  */
  myWatch->Stop();
  cout << myWatch->RealTime() << " seconds has passed! " << endl; 

}

void trcls(const char *fname, const char *oname, const char *tkdbc)
{
  AMSChain ch;
  TString sfn = fname;
  if (!sfn.Contains(".root")) sfn += "*.root";
  if (ch.Add(sfn) <= 0) return;

  Int_t ntr  = ch.GetNtrees();
  Int_t nent = ch.GetEntries();
  if (ntr <= 0 || nent <= 0) return;

  cout << "Ntr,Nent= " << ntr << " " << nent << endl;

  for (Int_t i = 0; i < ntr; i++)
    cout << ch.GetListOfFiles()->At(i)->GetTitle() << endl;

  Int_t   idata[13];
  Float_t fdata[28];

  TFile of(oname, "recreate");
  TTree *tree = new TTree("tree", "trcls");
  tree->Branch("idata", idata, "run/I:event/I:ient/I:time/I:tkml[9]/I");
  tree->Branch("fdata", fdata, "engc/F:enge/F:rgt/F:chrg/F:"
	                       "p0x/F:p0y/F:dzx/F:dzy/F:csqx/F:csqy/F:"
	                       "xcog[9]/F:ycog[9]/F");

  if (tkdbc && tkdbc[0] && tkdbc[0] != '0') {
    ch.GetEvent(0);
    TkDBc::Head->init(3, tkdbc);
  }
  TrExtAlignDB::OverLoadFlag = 0;

  Int_t malg = 2; // kAlcaraz | kMultScat
  Int_t pat0 = 3; // Inner only

  Int_t nevt = 0;
  Int_t nrsl = 0;
  Int_t npsl = 0;
  Int_t nfil = 0;
  Int_t intv = 10000;

  signal(SIGTERM, handler);
  signal(SIGINT,  handler);

  TStopwatch timer;
  timer.Start();

  for (Int_t ient = 0; ient < nent && !SigTERM; ient++) {
    AMSEventR *evt = ch.GetEvent(ient);
    nevt++;

    if (nevt%intv == 0 || nevt == nent) {
      Double_t tm = timer.RealTime();
      timer.Continue();
      cout << Form("%6d %6d %6d %7d (%5.1f%%) %4.0f sec (%4.1f kHz)",
		   nrsl, npsl, nfil, nevt,
		   100.*nevt/nent, tm, nevt/tm*1e-3)
	   << endl;
    }
    if (evt->nTrRecHit() >= 1600) continue;
    if (evt->nTrTrack () !=    1) continue;

    TrTrackR *trk = evt->pTrTrack(0);
    if (!trk) continue;

    Double_t chgp = TrCharge::GetQ(trk, 1);
    Double_t chgn = TrCharge::GetQ(trk, 0);
    if (chgp <= 0 || chgn <= 0) continue;

    Int_t mfp = trk->iTrTrackPar(malg, pat0, 0);
    if (mfp < 0) continue;

    Double_t rgtp = trk->GetRigidity(mfp);
    if (TMath::Abs(rgtp) < 5) continue;


    //////////////////// Recalc and refit ////////////////////
    trk->RecalcHitCoordinates();

    Int_t mf0 = trk->iTrTrackPar(malg, pat0, 2);
    if (mf0 < 0) continue;

    Double_t rgt0 = trk->GetRigidity(mf0);
    if (TMath::Abs(rgt0) < 10) continue;
    nrsl++;


    //////////////////// Pre-selection ////////////////////
    Bool_t psel = kTRUE;
    Int_t  span = (TrTrackSelection::GetSpanFlags(trk) & 0xff);
    if (!(span & TrTrackSelection::kMaxInt) ||
	!(span & TrTrackSelection::kAllPlane)) psel = kFALSE;

    if   (!(span & TrTrackSelection::kHalfL1N)) {
      if (!(span & TrTrackSelection::kHalfL9)) psel = kFALSE;

      AMSPoint pnt = trk->InterpolateLayerJ(9);
      if (TMath::Abs(pnt.x()) > 33) psel = kFALSE;
    }
    if (psel) npsl++;
//.........这里部分代码省略.........