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

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

		resultDir = TString(gSystem->Getenv("RESULT_DIR"));

		stsMatBudgetFileName = TString(gSystem->Getenv("STS_MATERIAL_BUDGET_FILE"));
		trdAnnCut = TString(gSystem->Getenv("TRD_ANN_CUT")).Atof();
		minNofPointsTrd = TString(gSystem->Getenv("MIN_NOF_POINTS_TRD")).Atof();
	}

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

   TStopwatch timer;
   timer.Start();

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

	// -----   Reconstruction run   -------------------------------------------
	FairRunAna *run= new FairRunAna();
	if (mcFile != "") run->SetInputFile(mcFile);
	if (recoFile != "") run->AddFriend(recoFile);
	if (qaFile != "") run->SetOutputFile(qaFile);

	CbmKF* kalman = new CbmKF();
	run->AddTask(kalman);
	CbmL1* l1 = new CbmL1();
	l1->SetMaterialBudgetFileName(stsMatBudgetFileName);
	run->AddTask(l1);

   // Reconstruction Qa
   CbmLitTrackingQa* trackingQa = new CbmLitTrackingQa();
   trackingQa->SetMinNofPointsSts(4);
   trackingQa->SetUseConsecutivePointsInSts(true);
   trackingQa->SetMinNofPointsTrd(minNofPointsTrd);
   trackingQa->SetMinNofPointsMuch(10);
   trackingQa->SetMinNofPointsTof(1);
   trackingQa->SetQuota(0.7);
   trackingQa->SetMinNofHitsTrd(minNofPointsTrd);
   trackingQa->SetMinNofHitsMuch(10);
   trackingQa->SetVerbose(0);
   trackingQa->SetMinNofHitsRich(7);
   trackingQa->SetQuotaRich(0.6);
   trackingQa->SetPRange(30, 0., 6.);
   trackingQa->SetOutputDir(std::string(resultDir));
   std::vector<std::string> trackCat, richCat;
   trackCat.push_back("All");
   trackCat.push_back("Electron");
   richCat.push_back("All");
   richCat.push_back("Electron");
   richCat.push_back("ElectronReference");
   trackingQa->SetTrackCategories(trackCat);
   trackingQa->SetRingCategories(richCat);
   trackingQa->SetTrdAnnCut(trdAnnCut);
   run->AddTask(trackingQa);

   CbmLitFitQa* fitQa = new CbmLitFitQa();
   fitQa->SetMvdMinNofHits(0);
   fitQa->SetStsMinNofHits(4);
   fitQa->SetMuchMinNofHits(10);
   fitQa->SetTrdMinNofHits(minNofPointsTrd);
   fitQa->SetPRange(30, 0., 3.);
   fitQa->SetOutputDir(std::string(resultDir));
   run->AddTask(fitQa);

  /* CbmLitClusteringQa* clusteringQa = new CbmLitClusteringQa();
   clusteringQa->SetMuchDigiFileName(muchDigiFile.Data());
   clusteringQa->SetOutputDir(std::string(resultDir));
   run->AddTask(clusteringQa);*/

   CbmLitTofQa* tofQa = new CbmLitTofQa();
   tofQa->SetOutputDir(std::string(resultDir));
   run->AddTask(tofQa);


    // -----  Parameter database   --------------------------------------------
   FairRuntimeDb* rtdb = run->GetRuntimeDb();
   FairParRootFileIo* parIo1 = new FairParRootFileIo();
   FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
   parIo1->open(parFile.Data());
   parIo2->open(parFileList, "in");
   rtdb->setFirstInput(parIo1);
   rtdb->setSecondInput(parIo2);
   rtdb->setOutput(parIo1);
   rtdb->saveOutput();

    run->Init();
    run->Run(0, nEvents);

    // -----   Finish   -------------------------------------------------------
    timer.Stop();
    std::cout << "Macro finished successfully." << std::endl;
    std::cout << "Output file is "    << recoFile << std::endl;
    std::cout << "Parameter file is " << parFile << std::endl;
    std::cout << "Real time " << timer.RealTime() << " s, CPU time " << timer.CpuTime() << " s" << std::endl;
    std::cout << " Test passed" << std::endl;
    std::cout << " All ok " << std::endl;
}

void run_trac_its(Int_t nEvents = 10, TString mcEngine = "TGeant3"){
        // Initialize logger
        FairLogger *logger = FairLogger::GetLogger();
        logger->SetLogVerbosityLevel("LOW");
        logger->SetLogScreenLevel("INFO");

        // Input and output file name
        std::stringstream inputfile, outputfile, paramfile;
        inputfile << "AliceO2_" << mcEngine << ".clus_" << nEvents << "_event.root";
        paramfile << "AliceO2_" << mcEngine << ".params_" << nEvents << ".root";
        outputfile << "AliceO2_" << mcEngine << ".trac_" << nEvents << "_event.root";

        // Setup timer
        TStopwatch timer;

        // Setup FairRoot analysis manager
        FairRunAna * fRun = new FairRunAna();
        FairFileSource *fFileSource = new FairFileSource(inputfile.str().c_str());
        fRun->SetSource(fFileSource);
        fRun->SetOutputFile(outputfile.str().c_str());

        // Setup Runtime DB
        FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
        FairParRootFileIo* parInput1 = new FairParRootFileIo();
        parInput1->open(paramfile.str().c_str());
        rtdb->setFirstInput(parInput1);

        // Setup tracker
        // To run with n threads call AliceO2::ITS::CookedTrackerTask(n)
        AliceO2::ITS::CookedTrackerTask *trac = new AliceO2::ITS::CookedTrackerTask;

        fRun->AddTask(trac);

        fRun->Init();

        AliceO2::Field::MagneticField* fld = (AliceO2::Field::MagneticField*)fRun->GetField();
      	if (!fld) {
      	  std::cout << "Failed to get field instance from FairRunAna" << std::endl;
      	  return;
      	}
      	trac->setBz(fld->solenoidField()); //in kG

        timer.Start();
        fRun->Run();

        std::cout << std::endl << std::endl;

        // Extract the maximal used memory an add is as Dart measurement
        // This line is filtered by CTest and the value send to CDash
        FairSystemInfo sysInfo;
        Float_t maxMemory=sysInfo.GetMaxMemory();
        std::cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">";
        std::cout << maxMemory;
        std::cout << "</DartMeasurement>" << std::endl;

        timer.Stop();
        Double_t rtime = timer.RealTime();
        Double_t ctime = timer.CpuTime();

        Float_t cpuUsage=ctime/rtime;
        cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
        cout << cpuUsage;
        cout << "</DartMeasurement>" << endl;
        cout << endl << endl;
        cout << "Macro finished succesfully." << endl;

        std::cout << endl << std::endl;
        std::cout << "Output file is "    << outputfile.str() << std::endl;
        //std::cout << "Parameter file is " << parFile << std::endl;
        std::cout << "Real time " << rtime << " s, CPU time " << ctime
                  << "s" << endl << endl;
}

void ana_Main_MC(TString ds="relval", TString physics="ttbar") {

	gSystem->Load("libSusyEvent.so");

	// Look ../jec/JetMETObjects/README
	gSystem->Load("../jec/lib/libJetMETObjects.so");

	// Printing utility for ntuple variables
	gROOT->LoadMacro("SusyEventPrinter.cc+");

	// Main analysis code
	gROOT->LoadMacro("SusyMainAna_MC.cc+");

	// chain of inputs
	TChain* chain = new TChain("susyTree");

	//////////////// MC files /////////////////
	cout<<"I survive this long1 "<< which_MC_to_use<< endl;
	MCpoint* thisMCpoint = setupMCpoint(which_MC_to_use);
	cout<<"I survive this long2"<<endl;
	chain->Add(thisMCpoint->filepath.c_str());
	cout<<"I survive this long"<<endl;


	//chain->Add("../susyEvents_AB_1M_ho200_v2.root");
	//chain->Add("../susyEvents_newNatural.root"); //last used!!
	//chain->Add("/eos/uscms/store/user/abarker/MC/newNat350_225/MC_AB_2500k_NEWnaturalHiggsinoNLSPout_mst_350_M3_5025_mu_225.root");//same thing as ../susyEvents_newNatural.root
	//chain->Add("/eos/uscms/store/user/abarker/MC/st_250_ho_150/MC_AB_2500k_st_250_ho_150.root");
	//chain->Add("/eos/uscms/store/user/abarker/MC/st_250_ho_200/MC_AB_2500k_st_250_ho_200.root");
	//chain->Add("/eos/uscms/store/user/abarker/MC/st_350_ho_200/MC_AB_2500k_mst_350_mu_200.root");
	//chain->Add("/eos/uscms/store/user/abarker/MC/ho_140/MC_AB_2500k_ho_140.root");
	//chain->Add("/eos/uscms/store/user/abarker/MC/ho_200/MC_AB_2500k_ho_200.root");


	//chain->Add("../susyEvents_newNatural.root");
	//chain->Add("dcache:/pnfs/cms/WAX/resilient/abarker/MC/MC_AB_2500k_NEWnaturalHiggsinoNLSPout_mst_350_M3_5025_mu_225.root");

	//chain->Add("dcache:/pnfs/cms/WAX/resilient/abarker/MC/MC_AB_2500k_st_250_ho_150.root");
	//chain->Add("dcache:/pnfs/cms/WAX/resilient/abarker/MC/MC_AB_2500k_st_250_ho_200.root");
	//chain->Add("dcache:/pnfs/cms/WAX/resilient/abarker/MC/MC_AB_2500k_mst_350_mu_200.root");
	//chain->Add("dcache:/pnfs/cms/WAX/resilient/abarker/MC/MC_AB_2500k_ho_140.root");
	//chain->Add("dcache:/pnfs/cms/WAX/resilient/abarker/MC/MC_AB_2500k_ho_200.root");


	SusyMainAna_MC* sea = new SusyMainAna_MC(chain);

	// configuration parameters
	// any values given here will replace the default values
	sea->SetDataset(physics+"_"+ds);        // dataset name
	sea->SetPrintInterval(1e4);             // print frequency
	sea->SetPrintLevel(0);                  // print level for event contents
	sea->SetUseTrigger(false);
	/*
	   sea->AddHltName("HLT_Photon36_CaloIdL_Photon22_CaloIdL");  // add HLT trigger path name 
	   sea->AddHltName("HLT_Photon32_CaloIdL_Photon26_CaloIdL");  // add HLT trigger path name
	   sea->AddHltName("HLT_Photon26_R9Id85_Photon18_R9Id85_Mass60");
	   sea->AddHltName("HLT_Photon26_R9Id85_Photon18_CaloId10_Iso50_Mass60");
	   sea->AddHltName("HLT_Photon26_CaloId10_Iso50_Photon18_R9Id85_Mass60");
	   sea->AddHltName("HLT_Photon26_CaloId10_Iso50_Photon18_CaloId10_Iso50_Mass60");
	   sea->AddHltName("HLT_Photon26_R9Id85_OR_CaloId10_Iso50_Photon18_R9Id85_OR_CaloId10_Iso50_Mass60");
	   sea->AddHltName("HLT_Photon26_R9Id85_OR_CaloId10_Iso50_Photon18_R9Id85_OR_CaloId10_Iso50_Mass70");
	   sea->AddHltName("HLT_Photon36_R9Id85_Photon22_R9Id85");
	   sea->AddHltName("HLT_Photon36_R9Id85_Photon22_CaloId10_Iso50");
	   sea->AddHltName("HLT_Photon36_CaloId10_Iso50_Photon22_R9Id85");
	   sea->AddHltName("HLT_Photon36_CaloId10_Iso50_Photon22_CaloId10_Iso50");
	   sea->AddHltName("HLT_Photon36_R9Id85_OR_CaloId10_Iso50_Photon22_R9Id85_OR_CaloId10_Iso50");
	 */
	sea->SetFilter(false);                  // filter events passing final cuts
	sea->SetProcessNEvents(-1);             // number of events to be processed

	// as an example -- add your favorite Json here.  More than one can be "Include"ed
	//  sea->IncludeAJson("Cert_161079-161352_7TeV_PromptReco_Collisions11_JSON_noESpbl_v2.txt");
	//sea->IncludeAJson("anotherJSON.txt");

	TStopwatch ts;

	ts.Start();

	sea->Loop();

	ts.Stop();

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

}

//.........这里部分代码省略.........
    // mTof
    run->AddModule(new R3BmTof("mtof_v17a.geo.root", { -155.824045, 0.523976, 761.870346 }, { "", -90., +16.7, 90. }));

    // MFI
    //run->AddModule(new R3BMfi("mfi_v17a.geo.root", { -63.82, 0., 520.25 }, { "", 90., +13.5, 90. })); // s412

    // NeuLAND
    // run->AddModule(new R3BNeuland("neuland_test.geo.root", { 0., 0., 1400. + 12 * 5. }));

    // -----   Create R3B  magnetic field ----------------------------------------
    // NB: <D.B>
    // If the Global Position of the Magnet is changed
    // the Field Map has to be transformed accordingly
    R3BGladFieldMap* magField = new R3BGladFieldMap("R3BGladMap");
    magField->SetScale(fieldScale);

    if (magnet == kTRUE)
    {
        run->SetField(magField);
    }
    else
    {
        run->SetField(NULL);
    }

    // -----   Create PrimaryGenerator   --------------------------------------
    // 1 - Create the Main API class for the Generator
    FairPrimaryGenerator* primGen = new FairPrimaryGenerator();

    if (generator.CompareTo("box") == 0)
    {
        FairIonGenerator* boxGen = new FairIonGenerator(50, 128, 50, 1, 0., 0., 1.3, 0., 0., 0.);
        primGen->AddGenerator(boxGen);
    }

    if (generator.CompareTo("ascii") == 0)
    {
        R3BAsciiGenerator* gen = new R3BAsciiGenerator((dir + "/input/" + inputFile).Data());
        primGen->AddGenerator(gen);
    }

    run->SetGenerator(primGen);

    run->SetStoreTraj(storeTrajectories);

    FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");
    FairLogger::GetLogger()->SetLogScreenLevel("INFO");

    // -----   Initialize simulation run   ------------------------------------
    run->Init();
    TVirtualMC::GetMC()->SetRandom(new TRandom3(randomSeed));

    // ------  Increase nb of step for CALO
    Int_t nSteps = -15000;
    TVirtualMC::GetMC()->SetMaxNStep(nSteps);

    // -----   Runtime database   ---------------------------------------------
    R3BFieldPar* fieldPar = (R3BFieldPar*)rtdb->getContainer("R3BFieldPar");
    if (NULL != magField)
    {
        fieldPar->SetParameters(magField);
        fieldPar->setChanged();
    }
    Bool_t kParameterMerged = kTRUE;
    FairParRootFileIo* parOut = new FairParRootFileIo(kParameterMerged);
    parOut->open(parFile.Data());
    rtdb->setOutput(parOut);
    rtdb->saveOutput();
    rtdb->print();

    // -----   Start run   ----------------------------------------------------
    if (nEvents > 0)
    {
        run->Run(nEvents);
    }

    // -----   Finish   -------------------------------------------------------
    timer.Stop();
    Double_t rtime = timer.RealTime();
    Double_t ctime = timer.CpuTime();
    cout << endl << endl;
    cout << "Macro finished succesfully." << endl;
    cout << "Output file is " << outFile << endl;
    cout << "Parameter file is " << parFile << endl;
    cout << "Real time " << rtime << " s, CPU time " << ctime << "s" << endl << endl;

    cout << " Test passed" << endl;
    cout << " All ok " << endl;

    // Snap a picture of the geometry
    // If this crashes, set "OpenGL.SavePicturesViaFBO: no" in your .rootrc
    /*gStyle->SetCanvasPreferGL(kTRUE);
    gGeoManager->GetTopVolume()->Draw("ogl");
    TGLViewer* v = (TGLViewer*)gPad->GetViewer3D();
    v->SetStyle(TGLRnrCtx::kOutline);
    v->RequestDraw();
    v->SavePicture("run_sim-side.png");
    v->SetPerspectiveCamera(TGLViewer::kCameraPerspXOZ, 25., 0, 0, -90. * TMath::DegToRad(), 0. * TMath::DegToRad());
    v->SavePicture("run_sim-top.png");*/
}

void Check_PFelecFix(int radius = 3,
		     char * algo = (char*)"PF",
		     char * bkgsub = (char*)"Pu"){

  
  TH1::SetDefaultSumw2();
  //gStyle->SetOptStat(0);
  
  TStopwatch timer;
  timer.Start();
  
  TDatime date;

  char * fileType[4][256] = {"badFile","Fix_1","Fix_2","Fix_3"};
  
  // get the input hiForest files, these are going to be in an array 0 - bad fine, 1 - fix_1, 2 - fix_2, 3 - fix_3; 
  TFile * fIn[4];

  // input the filen names
  fIn[0] = TFile::Open("");
  fIn[1] = TFile::Open("");
  fIn[2] = TFile::Open("");
  fIn[3] = TFile::Open("");
  
  // get the jet trees from the necessary files (these events already passed the event quality cuts so only hav to get the jet Tree).
  TTree * jet[4];
   
  // get the histograms from these files for the different centrality bins and the fixes. 
  TH1F * heMax[4][nbins_cent], * heSum[4][nbins_cent];
  TH2F * heMax_vs_jtpt[4][nbins_cent], * heMaxJtpt_vs_jtpt[4][nbins_cent], * heMaxSumcand_vs_jtpt[4][nbins_cent];


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

    jet[k] = (TTree*)fIn[k]->Get(Form("ak%s%d%sJetAnalyzer/t",bkgsub,radius,algo));

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

      heMax[k][i] = new TH1F(Form("heMax_%s_cent%d",fileType,i),"",200,0,200);
      heMax_vs_jtpt[k][i] = new TH1F(Form("heMax_vs_jtpt_%s_cent0",fileType,i),"",400,0,400,100,0,200);
      heMaxJtpt_vs_jtpt[k][i] = new TH1F(Form("heMaxJtpt_vs_jtpt_%s_cent0",fileType,i),"",400,0,400,100,0,10);
      heMaxSumcand_vs_jtpt[k][i] = new TH1F(Form("heMaxSumcand_vs_jtpt_%s_cent0",fileType,i),"",400,0,400,100,0,10);

      jet[k]->Draw(Form("eMax>>heMax_%s_cent%d",fileType,i),Form("hiBin>=%d && hiBin<%d",5 * boundaries_cent[i], 5 * boundaries_cent[i+1]),"goff");
      jet[k]->Draw(Form("eMax:jtpt>>heMax_vs_jtpt_%s_cent%d",fileType,i),Form("hiBin>=%d && hiBin<%d",5 * boundaries_cent[i], 5 * boundaries_cent[i+1]),"goff");
      jet[k]->Draw(Form("eMax/jtpt:jtpt>>heMaxJtpt_vs_jtpt_%s_cent%d",fileType,i),Form("hiBin>=%d && hiBin<%d",5 * boundaries_cent[i], 5 * boundaries_cent[i+1]),"goff");
      jet[k]->Draw(Form("eMax/(chSum+neSum+muSum+phSum):jtpt>>heMaxSumcand_vs_jtpt_%s_cent%d",fileType,i),Form("hiBin>=%d && hiBin<%d",5 * boundaries_cent[i], 5 * boundaries_cent[i+1]),"goff");
    
    }

  }

  // now that we have the histograms, lets do the plotting part here.
  // we need one canvas for each plotted variable.
  // at the moment im only going to plot the most central events. 0 <= hiBin < 5

  TCanvas * ceMax, *ceMax_vs_jtpt, *ceMaxJtpt_vs_jtpt, *ceMaxSumcand_vs_jtpt;

  ceMax = new TCanvas("ceMax","",1200,1000);
  ceMax->Divide(4,1);
  ceMax_vs_jtpt = new TCanvas("ceMax_vs_jtpt","",1200,1000);
  ceMax_vs_jtpt->Divide(4,1);
  ceMaxJtpt_vs_jtpt = new TCanvas("ceMaxJtpt_vs_jtpt","",1200,1000);
  ceMaxJtpt_vs_jtpt->Divide(4,1);
  ceMaxSumcand_vs_jtpt = new TCanvas("ceMaxSumcand_vs_jtpt","",1200,1000);
  ceMaxSumcand_vs_jtpt->Divide(4,1);

  int cent = 0; // change this to draw other centrality classes 
  
  for(int k = 0; k<4; ++k){

    ceMax->cd(k+1);
    heMax[k][cent]->Draw();
    drawText(Form("%s %2.0f-%2.0f%",fileType,2.5*boundaries_cent[i],2.5*boundaries_cent[i+1]),0.2,0.2,14);

    ceMax_vs_jtpt->cd(k+1);
    heMax_vs_jtpt[k][cent]->Draw("colz");
    drawText(Form("%s %2.0f-%2.0f%",fileType,2.5*boundaries_cent[i],2.5*boundaries_cent[i+1]),0.2,0.2,14);

    ceMaxJtpt_vs_jtpt->cd(k+1);
    heMaxJtpt_vs_jtpt[k][cent]->Draw("colz");
    drawText(Form("%s %2.0f-%2.0f%",fileType,2.5*boundaries_cent[i],2.5*boundaries_cent[i+1]),0.2,0.2,14);

    ceMaxSumcand_vs_jtpt->cd(k+1);
    heMax_vs_jtpt[k][cent]->Draw("colz");
    drawText(Form("%s %2.0f-%2.0f%",fileType,2.5*boundaries_cent[i],2.5*boundaries_cent[i+1]),0.2,0.2,14);

  }

  ceMax->SaveAs(Form("PbPb_eMaxvariable_cent%d_ak%s%d%s_%d.pdf",cent,bkgsub,radius,algo,date.GetDate()),"RECREATE");
  ceMax_vs_jtpt->SaveAs(Form("PbPb_eMax_vs_jtpt_cent%d_ak%s%d%s_%d.pdf",cent,bkgsub,radius,algo,date.GetDate()),"RECREATE");
  ceMaxJtpt_vs_jtpt->SaveAs(Form("PbPb_eMaxOverJtpt_vs_jtpt_cent%d_ak%s%d%s_%d.pdf",cent,bkgsub,radius,algo,date.GetDate()),"RECREATE");
  ceMaxSumcand_vs_jtpt->SaveAs(Form("PbPb_eMaxOver_SumCandidates_without_eMax_vs_jtpt_cent%d_ak%s%d%s_%d.pdf",cent,bkgsub,radius,algo,date.GetDate()),"RECREATE");
  
  timer.Stop();
  cout<<"Macro finished: "<<endl;
  cout<<"CPU time (min)  = "<<(float)timer.CpuTime()/60<<endl;
  cout<<"Real time (min) = "<<(float)timer.RealTime()/60<<endl;


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

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 CrapPhoAnalysis(TString outputName="analysis"){

  // REMOVE THE LINE BELOW IF NOT RUNNING IN CMSSW ENVIRONMENT
  //gSystem->Load("libCondFormatsJetMETObjects.so");

  //gSystem->Load("../../SUSYPhotonAnalysis/SusyNtuplizer/macro/libSusyEvent.so");

  //gSystem->AddIncludePath("-I" + TString(gSystem->Getenv("CMSSW_RELEASE_BASE")) + "/src");

  // Analysis macro
  //gROOT->LoadMacro("SusyEventAnalyzer.cc+");

  //gSystem->AddIncludePath("-I /Users/dmason/play/photons/SUSYPhotonAnalysis/SusyNtuplizer/src");
  //gSystem->AddIncludePath("-I /Users/dmason/play/photons/Crap");

  //gSystem->Load("/Users/dmason/play/photons/SUSYPhotonAnalysis/SusyNtuplizer/src/SusyEvent.h+");
  gROOT->ProcessLine(".L ../SusyEvent.h+");
  gROOT->ProcessLine(".L ../Crap.h+");
  
  gROOT->ProcessLine(".L CrapPhoAnalysis.h+");


  gSystem->Load("../lib/libCrapPhoAnalysis.so");
  
  //gROOT->LoadMacro("CrapPhoAnalysis.cc++");

  // chain of inputs
  TChain chain("susyTree");
  //chain.Add("susyEvents.root");
  chain.Add("root://cmseos.fnal.gov//eos/uscms/store/user/lpcpjm/SusyNtuples/cms538v1/Run2012D-22Jan2013-v1/DoublePhoton/susyEvents_361_1_fC3.root");

  // Disabling unused branches will speed up the processing significantly, but risks inconsistencies if wrong trees are turned off.
  // Make sure you know what you are doing.
  // Especially be careful when producing a skim - disabled branches will not be copied.
  // You can either enable all branches, or have two Event objects (one for event selection with reduced branch configuration, and the
  // other for event copying with branches fully enabled).
  //chain.SetBranchStatus("*", 0);
  //chain.SetBranchStatus("runNumber", 1);
  //chain.SetBranchStatus("luminosityBlockNumber", 1);
  //chain.SetBranchStatus("eventNumber", 1);
  //chain.SetBranchStatus("isRealData", 1);
  //chain.SetBranchStatus("metFilterBit", 1);
  //chain.SetBranchStatus("rho", 1);
  //chain.SetBranchStatus("rho25", 1);
  //chain.SetBranchStatus("hlt*", 1);
  //chain.SetBranchStatus("vertices*", 1);
  //chain.SetBranchStatus("photons_photons*", 1);
  //chain.SetBranchStatus("muons_muons*", 1);
  //chain.SetBranchStatus("electrons_gsfElectrons*", 1);
  //chain.SetBranchStatus("pfJets_ak5*", 1);
  //chain.SetBranchStatus("met_pfType01CorrectedMet*", 1);
  //chain.SetBranchStatus("pfParticles*", 1);

  if(chain.LoadTree(0) != 0){
    cerr << "Error with input chain. Do the files exist?" << endl;
    return;
  }

  CrapPhoAnalysis sea(chain);

  sea.SetOutput(outputName);
  sea.SetLogFile("cout"); // set to a full path to a file to output log to a file
  sea.SetPrintInterval(1000);
  sea.SetPrintLevel(0);
  sea.AddHltName("HLT_Photon36_CaloId10_Iso50_Photon22_CaloId10_Iso50");
  sea.CopyEvents(true);
  sea.SetProcessNEvents(-1);

  TStopwatch ts;

  ts.Start();

  sea.Run();

  ts.Stop();

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

}

//.........这里部分代码省略.........
	// add the gamma generator
	primGen->AddGenerator(gammasGen);
  } 
	

 if (fGenerator.CompareTo("r3b") == 0  ) {
  R3BSpecificGenerator *pR3bGen = new R3BSpecificGenerator();

  // R3bGen properties
  pR3bGen->SetBeamInteractionFlag("off");
  pR3bGen->SetBeamInteractionFlag("off");
  pR3bGen->SetRndmFlag("off");
  pR3bGen->SetRndmEneFlag("off");
  pR3bGen->SetBoostFlag("off");
  pR3bGen->SetReactionFlag("on");
  pR3bGen->SetGammasFlag("off");
  pR3bGen->SetDecaySchemeFlag("off");
  pR3bGen->SetDissociationFlag("off");
  pR3bGen->SetBackTrackingFlag("off");
  pR3bGen->SetSimEmittanceFlag("off");

  // R3bGen Parameters
  pR3bGen->SetBeamEnergy(1.); // Beam Energy in GeV
  pR3bGen->SetSigmaBeamEnergy(1.e-03); // Sigma(Ebeam) GeV
  pR3bGen->SetParticleDefinition(2212); // Use Particle Pdg Code
  pR3bGen->SetEnergyPrim(0.3); // Particle Energy in MeV
  Int_t fMultiplicity = 50;
  pR3bGen->SetNumberOfParticles(fMultiplicity); // Mult.

  // Reaction type
  //        1: "Elas"
  //        2: "iso"
  //        3: "Trans"
  pR3bGen->SetReactionType("Elas");

  // Target  type
  //        1: "LeadTarget"
  //        2: "Parafin0Deg"
  //        3: "Parafin45Deg"
  //        4: "LiH"

  pR3bGen->SetTargetType(Target.Data());
  Double_t thickness = (0.11/2.)/10.;  // cm
  pR3bGen->SetTargetHalfThicknessPara(thickness); // cm
  pR3bGen->SetTargetThicknessLiH(3.5); // cm
  pR3bGen->SetTargetRadius(1.); // cm

  pR3bGen->SetSigmaXInEmittance(1.); //cm
  pR3bGen->SetSigmaXPrimeInEmittance(0.0001); //cm

  // Dump the User settings
  pR3bGen->PrintParameters();  
  primGen->AddGenerator(pR3bGen);
  }


  run->SetGenerator(primGen);


  //-------Set visualisation flag to true------------------------------------
  if (fVis==kTRUE){
     run->SetStoreTraj(kTRUE);
  }else{
     run->SetStoreTraj(kFALSE);
  }   

  // -----   Initialize simulation run   ------------------------------------
  run->Init();

  // ------  Increase nb of step for CALO
  Int_t nSteps = -15000;
  gMC->SetMaxNStep(nSteps);

  // -----   Runtime database   ---------------------------------------------
  Bool_t kParameterMerged = kTRUE;
  FairParRootFileIo* parOut = new FairParRootFileIo(kParameterMerged);
  parOut->open(ParFile.Data());
  rtdb->setOutput(parOut);
  rtdb->saveOutput();
  rtdb->print();
   
  // -----   Start run   ----------------------------------------------------
  if (nEvents>0) run->Run(nEvents);
  
  // -----   Finish   -------------------------------------------------------
  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
  cout << endl << endl;
  cout << "Macro finished succesfully." << endl;
  cout << "Output file is "    << OutFile << endl;
  cout << "Parameter file is " << ParFile << endl;
  cout << "Real time " << rtime << " s, CPU time " << ctime 
       << "s" << endl << endl;
  // ------------------------------------------------------------------------

  cout << " Test passed" << endl;
  cout << " All ok " << endl;

}

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


                int ieta=-1;
                if(fabs(recoeta)<1.3)ieta=0; //! barrel region
                else ieta=1; //! HCAL region


                //! Response in eta
                hratiocorrrefpt_eta[nj][icent][ieta]->Fill(refpt[gj],recopt/refpt[gj],wxs*wcen*wvz);
                hratiorawrefpt_eta [nj][icent][ieta]->Fill(refpt[gj],rawpt[gj]/refpt[gj],wxs*wcen*wvz);

                //! Ratio of recopT / rawpT
                hrajptrpt[nj][icent] ->Fill(recopt,recopt/rawpt[gj],wxs*wcen*wvz);

                //! Jet eta, phi, pt, eta-pt, eta-phi and pt-phi
                hjeteta   [nj][icent]->Fill(jteta[gj],wxs*wcen*wvz);
                hjetphi   [nj][icent]->Fill(jtphi[gj],wxs*wcen*wvz);
                hjetpteta [nj][icent]->Fill(recopt,jteta[gj],wxs*wcen*wvz);
                hjetptphi [nj][icent]->Fill(jtpt[gj],jtphi[gj],wxs*wcen*wvz);
                hjetetaphi[nj][icent]->Fill(jteta[gj],jtphi[gj],wxs*wcen*wvz);

                hgenpt_genm [nj][icent]->Fill(refpt[gj],wxs*wcen*wvz);
                hrecopt_genm[nj][icent]->Fill(recopt,wxs*wcen*wvz);
                hrawpt_genm [nj][icent]->Fill(rawpt[gj],wxs*wcen*wvz);

                //! Very fine bin in ref pt used for response
                hrescrpt_genm [nj][icent]->Fill(refpt[gj],recopt/refpt[gj],wxs*wcen*wvz);
                hresrrpt_genm [nj][icent]->Fill(refpt[gj],rawpt[gj]/refpt[gj],wxs*wcen*wvz);
                hresrcrpt_genm[nj][icent]->Fill(recopt,recopt/rawpt[gj],wxs*wcen*wvz);

                //! Very fine bin in ref eta used for response
                hrescreta_genm [nj][icent]->Fill(refeta[gj],refpt[gj],recopt/refpt[gj],wxs*wcen*wvz);
                hresrreta_genm [nj][icent]->Fill(refeta[gj],refpt[gj],rawpt[gj]/refpt[gj],wxs*wcen*wvz);
                hresrcreta_genm[nj][icent]->Fill(recoeta,refpt[gj],recopt/rawpt[gj],wxs*wcen*wvz);

                //! Very fine bin in ref phi used for response
                hrescrphi_genm [nj][icent]->Fill(refphi[gj],refpt[gj],recopt/refpt[gj],wxs*wcen*wvz);
                hresrrphi_genm [nj][icent]->Fill(refphi[gj],refpt[gj],rawpt[gj]/refpt[gj],wxs*wcen*wvz);
                hresrcrphi_genm[nj][icent]->Fill(recophi,refpt[gj],recopt/rawpt[gj],wxs*wcen*wvz);

                //! Response in different eta and phi bins
                int etabin = GetEtaBin(fabs(refeta[gj]));
                int phibin = GetPhiBin(refphi[gj]);

                //! Response in eta and phi bins
                if(etabin >= 0 && etabin<maxe) {
                    hpteta[nj][icent][etabin]->Fill(refpt[gj],recopt/refpt[gj],wxs*wcen*wvz);
                }
                if(phibin >= 0 && phibin<maxph) {
                    hptphi[nj][icent][phibin]->Fill(refpt[gj],recopt/refpt[gj],wxs*wcen*wvz);
                }

            }//! igen loop


            iEvent++;
            //std::cout<<"Completed event #  "<<ievt<<std::endl;

        }//! event loop ends

        //delete parHI_l2;
        //delete parHI_l3;
        //delete _JEC_HI;

    }//! jet loop

    std::cout<<std::endl;
    for(int nj=0; nj<knj; nj++) {
        std::cout<<calgo[nj] << std::endl;
        if(strcmp(ksp,"pp")==0) {
            std::cout<<"\t # of events : "<< " pp " << "  " << nevt[nj][0] << "  " << njets[nj][0] << std::endl;
        }
        else {
            for(int ic=0; ic<ncent; ic++) {
                std::cout<<"\t # of events : "<< ccent[ic] << "  " << nevt[nj][ic] << "  " << njets[nj][ic] << std::endl;
            }
        }
        std::cout<<std::endl;
    }
    std::cout<<std::endl;



    //! Write to output file
    fout->cd();
    fout->Write();
    fout->Close();


    //! Check
    timer.Stop();
    double rtime  = timer.RealTime();
    double ctime  = timer.CpuTime();

    std::cout<<std::endl;
    std::cout<<Form("RealTime=%f seconds, CpuTime=%f seconds",rtime,ctime)<<std::endl;
    std::cout<<std::endl;
    std::cout<<"Good bye : " <<"\t"<<std::endl;
    return 1;
}

//.........这里部分代码省略.........
    relrespcorr2->SetMarkerColor(3);
    relrespcorr2->SetLineColor(3);
    relrespcorr3->SetMarkerStyle(29);
    relrespcorr3->SetMarkerColor(4);
    relrespcorr3->SetLineColor(4);
    relrespcorr4->SetMarkerStyle(33);
    relrespcorr4->SetMarkerColor(1);
    relrespcorr4->SetLineColor(1);
    relrespcorr5->SetMarkerStyle(34);
    relrespcorr5->SetMarkerColor(7);
    relrespcorr5->SetLineColor(7);
    relrespcorr1->Draw("AP");
    relrespcorr2->Draw("P");
    relrespcorr3->Draw("P");
    relrespcorr4->Draw("P");
    relrespcorr5->Draw("P");
    TLegend *legendc = new TLegend(0.4,0.6,0.89,0.89);
    legendc->AddEntry((TObject*)0, "Corrected pp Data Relative Response", "");
    legendc->AddEntry((TObject*)0, Form("ak%iPFJetAnalyzer",radius), "");
    legendc->AddEntry(relrespcorr5, "40 < p^{avg}_{T} < 60 GeV/c", "p");
    legendc->AddEntry(relrespcorr4, "60 < p^{avg}_{T} < 80 GeV/c", "p");
    legendc->AddEntry(relrespcorr2, "80 < p^{avg}_{T} < 100 GeV/c", "p");
    legendc->AddEntry(relrespcorr1, "100 < p^{avg}_{T} < 140 GeV/c", "p");
    legendc->AddEntry(relrespcorr3, "140 < p^{avg}_{T} < 200 GeV/c", "p");
    legendc->AddEntry((TObject*)0, "|#Delta#phi| > 2.5", "");
    legendc->AddEntry((TObject*)0, "#alpha = p^{third}_{T}/p^{avg}_{T} < 0.2", "");
    legendc->AddEntry((TObject*)0, "|#eta^{ref}| < 1.3", "");
    legendc->AddEntry((TObject*)0, "|#eta^{probe}| < 3.0", "");
    legendc->Draw();
    c->Update();

    TCanvas *d = new TCanvas();
    casymcorr1->SetMarkerStyle(21);
    casymcorr1->SetMarkerColor(2);
    casymcorr1->SetLineColor(2);
    casymcorr2->SetMarkerStyle(22);
    casymcorr2->SetMarkerColor(3);
    casymcorr2->SetLineColor(3);
    casymcorr3->SetMarkerStyle(29);
    casymcorr3->SetMarkerColor(4);
    casymcorr3->SetLineColor(4);
    casymcorr4->SetMarkerStyle(33);
    casymcorr4->SetMarkerColor(1);
    casymcorr4->SetLineColor(1);
    casymcorr5->SetMarkerStyle(34);
    casymcorr5->SetMarkerColor(7);
    casymcorr5->SetLineColor(7);
    casymcorr1->Draw("AP");
    casymcorr2->Draw("P");
    casymcorr3->Draw("P");
    casymcorr4->Draw("P");
    casymcorr5->Draw("P");
    TLegend *legendd = new TLegend(0.4,0.6,0.89,0.89);
    legendd->AddEntry((TObject*)0, "Asymmetry Correction Factor", "");
    legendd->AddEntry((TObject*)0, Form("ak%iPFJetAnalyzer",radius), "");
    legendd->AddEntry(casymcorr5, "40 < p^{avg}_{T} < 60 GeV/c", "p");
    legendd->AddEntry(casymcorr4, "60 < p^{avg}_{T} < 80 GeV/c", "p");
    legendd->AddEntry(casymcorr2, "80 < p^{avg}_{T} < 100 GeV/c", "p");
    legendd->AddEntry(casymcorr1, "100 < p^{avg}_{T} < 140 GeV/c", "p");
    legendd->AddEntry(casymcorr3, "140 < p^{avg}_{T} < 200 GeV/c", "p");
    legendd->AddEntry((TObject*)0, "|#Delta#phi| > 2.5", "");
    legendd->AddEntry((TObject*)0, "#alpha = p^{third}_{T}/p^{avg}_{T} < 0.2", "");
    legendd->AddEntry((TObject*)0, "|#eta^{ref}| < 1.3", "");
    legendd->AddEntry((TObject*)0, "|#eta^{probe}| < 3.0", "");
    legendd->Draw();
    d->Update();


    myfile->cd();
    a->Write();
    b->Write();
    c->Write();
    d->Write();
    myfile->Close();
    montecarlo1->Close();
    montecarlo2->Close();
    montecarlo3->Close();
    montecarlo4->Close();
    montecarlo5->Close();
    ppdata1->Close();
    ppdata2->Close();
    ppdata3->Close();
    ppdata4->Close();
    ppdata5->Close();
    corr1->Close();
    corr2->Close();
    corr3->Close();
    corr4->Close();
    corr5->Close();
    casym1->Close();
    casym2->Close();
    casym3->Close();
    casym4->Close();
    casym5->Close();
  }//END RADIUS LOOP
  timer.Stop();
  cout << "End of Macro Reached" << endl;
  cout << "CPU Time (sec)  : " << timer.CpuTime() << endl;
  cout << "Real Time (sec) : " << timer.RealTime() << endl;
}

//.........这里部分代码省略.........
chain->Add("looseSkims/SuperSkim_data_12C4_10.root");
chain->Add("looseSkims/SuperSkim_data_12C4_11.root");
chain->Add("looseSkims/SuperSkim_data_12C4_12.root");
chain->Add("looseSkims/SuperSkim_data_12C4_2.root");
chain->Add("looseSkims/SuperSkim_data_12C4_3.root");
chain->Add("looseSkims/SuperSkim_data_12C4_4.root");
chain->Add("looseSkims/SuperSkim_data_12C4_5.root");
chain->Add("looseSkims/SuperSkim_data_12C4_6.root");
chain->Add("looseSkims/SuperSkim_data_12C4_7.root");
chain->Add("looseSkims/SuperSkim_data_12C4_8.root");
chain->Add("looseSkims/SuperSkim_data_12C4_9.root");
}


if(useC5){
chain->Add("looseSkims/SuperSkim_data_12C5_0.root");
chain->Add("looseSkims/SuperSkim_data_12C5_1.root");
chain->Add("looseSkims/SuperSkim_data_12C5_10.root");
chain->Add("looseSkims/SuperSkim_data_12C5_11.root");
chain->Add("looseSkims/SuperSkim_data_12C5_12.root");
chain->Add("looseSkims/SuperSkim_data_12C5_13.root");
chain->Add("looseSkims/SuperSkim_data_12C5_14.root");
chain->Add("looseSkims/SuperSkim_data_12C5_15.root");
chain->Add("looseSkims/SuperSkim_data_12C5_16.root");
chain->Add("looseSkims/SuperSkim_data_12C5_17.root");
chain->Add("looseSkims/SuperSkim_data_12C5_18.root");
chain->Add("looseSkims/SuperSkim_data_12C5_19.root");
chain->Add("looseSkims/SuperSkim_data_12C5_2.root");
chain->Add("looseSkims/SuperSkim_data_12C5_20.root");
chain->Add("looseSkims/SuperSkim_data_12C5_21.root");
chain->Add("looseSkims/SuperSkim_data_12C5_22.root");
chain->Add("looseSkims/SuperSkim_data_12C5_23.root");
chain->Add("looseSkims/SuperSkim_data_12C5_24.root");
chain->Add("looseSkims/SuperSkim_data_12C5_25.root");
chain->Add("looseSkims/SuperSkim_data_12C5_26.root");
chain->Add("looseSkims/SuperSkim_data_12C5_27.root");
chain->Add("looseSkims/SuperSkim_data_12C5_28.root");
chain->Add("looseSkims/SuperSkim_data_12C5_29.root");
chain->Add("looseSkims/SuperSkim_data_12C5_3.root");
chain->Add("looseSkims/SuperSkim_data_12C5_30.root");
chain->Add("looseSkims/SuperSkim_data_12C5_31.root");
chain->Add("looseSkims/SuperSkim_data_12C5_32.root");
chain->Add("looseSkims/SuperSkim_data_12C5_33.root");
chain->Add("looseSkims/SuperSkim_data_12C5_4.root");
chain->Add("looseSkims/SuperSkim_data_12C5_5.root");
chain->Add("looseSkims/SuperSkim_data_12C5_6.root");
chain->Add("looseSkims/SuperSkim_data_12C5_7.root");
chain->Add("looseSkims/SuperSkim_data_12C5_8.root");
chain->Add("looseSkims/SuperSkim_data_12C5_9.root");
}
//if(useb1_dcache){
//chain->Add("dcache:/pnfs/cms/WAX/resilient/abarker/SusyLooseSkim_v2_2012/SuperSkim_data_12B1_37.root");
//chain->Add("dcache:/pnfs/cms/WAX/resilient/abarker/SusyLooseSkim_v2_2012/SuperSkim_data_12B1_79.root");
//}

  SusySkimmer* sea = new SusySkimmer(chain);

  // configuration parameters
  // any values given here will replace the default values
  sea->SetDataset(physics+"_"+ds);        // dataset name
  sea->SetPrintInterval(1e4);             // print frequency
  sea->SetPrintLevel(0);                  // print level for event contents
  sea->SetUseTrigger(true);
/*
  sea->AddHltName("HLT_Photon36_CaloIdL_Photon22_CaloIdL");  // add HLT trigger path name 
  sea->AddHltName("HLT_Photon32_CaloIdL_Photon26_CaloIdL");  // add HLT trigger path name
  sea->AddHltName("HLT_Photon26_R9Id85_Photon18_R9Id85_Mass60");
  sea->AddHltName("HLT_Photon26_R9Id85_Photon18_CaloId10_Iso50_Mass60");
  sea->AddHltName("HLT_Photon26_CaloId10_Iso50_Photon18_R9Id85_Mass60");
  sea->AddHltName("HLT_Photon26_CaloId10_Iso50_Photon18_CaloId10_Iso50_Mass60");
  sea->AddHltName("HLT_Photon26_R9Id85_OR_CaloId10_Iso50_Photon18_R9Id85_OR_CaloId10_Iso50_Mass60");
  sea->AddHltName("HLT_Photon26_R9Id85_OR_CaloId10_Iso50_Photon18_R9Id85_OR_CaloId10_Iso50_Mass70");
  sea->AddHltName("HLT_Photon36_R9Id85_Photon22_R9Id85");
  sea->AddHltName("HLT_Photon36_R9Id85_Photon22_CaloId10_Iso50");
  sea->AddHltName("HLT_Photon36_CaloId10_Iso50_Photon22_R9Id85");
  sea->AddHltName("HLT_Photon36_CaloId10_Iso50_Photon22_CaloId10_Iso50");
  sea->AddHltName("HLT_Photon36_R9Id85_OR_CaloId10_Iso50_Photon22_R9Id85_OR_CaloId10_Iso50");
*/
  sea->SetFilter(false);                  // filter events passing final cuts
  sea->SetProcessNEvents(-1);             // number of events to be processed

  // as an example -- add your favorite Json here.  More than one can be "Include"ed
  //sea->IncludeAJson("Cert_190456-199429_8TeV_PromptReco_Collisions12_JSON.txt");
  //sea->IncludeAJson("Cert_190456-202016_8TeV_PromptReco_Collisions12_JSON.txt");
  //sea->IncludeAJson("Cert_190456-203853_8TeV_PromptReco_Collisions12_JSON.txt");
  sea->IncludeAJson("Cert_190456-208686_8TeV_Collisions12_Abarker_Combined_JSON.txt");//combined 2012
  //sea->IncludeAJson("anotherJSON.txt");

  TStopwatch ts;

  ts.Start();

  sea->Loop();

  ts.Stop();

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

}

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

  TString outparam = "data/testparams_";
  outparam = outparam + mcEngine + ".root";

  
  fRun->SetOutputFile(outfile);


  // -----   Magnetic field   -------------------------------------------
  // Constant Field
  FairConstField *fMagField = new FairConstField();
  fMagField->SetField(0., 10. ,0. ); // values are in kG
  fMagField->SetFieldRegion(-50, 50,-50, 50, 350, 450);// values are in cm (xmin,xmax,ymin,ymax,zmin,zmax)
  fRun->SetField(fMagField);
  // --------------------------------------------------------------------

 
  // Set Material file Name
  //-----------------------
  fRun->SetMaterials("media.geo");
  
  // Create and add detectors
  //-------------------------
  FairModule *Cave= new FairCave("CAVE");
  Cave->SetGeometryFileName("cave.geo");
  fRun->AddModule(Cave);

  FairModule *Magnet= new FairMagnet("MAGNET");
  Magnet->SetGeometryFileName("magnet.geo");
  fRun->AddModule(Magnet);

  FairDetector *Torino= new FairTestDetector("TORINO", kTRUE);
  Torino->SetGeometryFileName("torino.geo");
  fRun->AddModule(Torino);

  // Create and Set Event Generator
  //-------------------------------

  FairPrimaryGenerator* primGen = new FairPrimaryGenerator();
  fRun->SetGenerator(primGen);
  
  // Box Generator
  FairBoxGenerator* boxGen = new FairBoxGenerator(2212, 10); // 13 = muon; 1 = multipl.
  boxGen->SetPRange(2., 2.); // GeV/c //setPRange vs setPtRange
  boxGen->SetPhiRange(0, 360); // Azimuth angle range [degree]
  boxGen->SetThetaRange(3, 10); // Polar angle in lab system range [degree]
  boxGen->SetCosTheta();//uniform generation on all the solid angle(default)

  // boxGen->SetXYZ(0., 0.37, 0.); 
  primGen->AddGenerator(boxGen);
  
  fRun->SetStoreTraj(kTRUE);
  
  fRun->Init();


 // -Trajectories Visualization (TGeoManager Only )
 // -----------------------------------------------

 // Set cuts for storing the trajectpries
  /* FairTrajFilter* trajFilter = FairTrajFilter::Instance();
     trajFilter->SetStepSizeCut(0.01); // 1 cm
     trajFilter->SetVertexCut(-2000., -2000., 4., 2000., 2000., 100.);
     trajFilter->SetMomentumCutP(10e-3); // p_lab > 10 MeV
     trajFilter->SetEnergyCut(0., 1.02); // 0 < Etot < 1.04 GeV
     trajFilter->SetStorePrimaries(kTRUE);
     trajFilter->SetStoreSecondaries(kTRUE);
   */


  // Fill the Parameter containers for this run
  //-------------------------------------------

  FairRuntimeDb *rtdb=fRun->GetRuntimeDb();
  Bool_t kParameterMerged=kTRUE;
  FairParRootFileIo* output=new FairParRootFileIo(kParameterMerged);
  output->open(outparam);
  rtdb->setOutput(output);

  rtdb->saveOutput();
  rtdb->print();

  // Transport nEvents
  // -----------------

  //  Int_t nEvents = 1;
  fRun->Run(nEvents);

  TString geoFile = tutdir + "/data/geofile_full.root";  
  fRun->CreateGeometryFile(geoFile.Data());


  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
  printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);

  cout << "Macro finished successfully." << endl;  

}  

void run_sim(Int_t nEvents = 100, TString mcEngine = "TGeant4")
{
    
  // Output file name
  TString outFile ="test.root";
    
  // Parameter file name
  TString parFile="params.root";
  
  // -----   Timer   --------------------------------------------------------
  TStopwatch timer;
  timer.Start();
  // ------------------------------------------------------------------------

  // -----   Create simulation run   ----------------------------------------
  FairRunSim* run = new FairRunSim();
  run->SetName(mcEngine);              // Transport engine
  run->SetOutputFile(new FairRootFileSink(outFile));          // Output file
  FairRuntimeDb* rtdb = run->GetRuntimeDb();
  // ------------------------------------------------------------------------
  
  // -----   Create media   -------------------------------------------------
  run->SetMaterials("media.geo");       // Materials
  // ------------------------------------------------------------------------
  
  // -----   Create geometry   ----------------------------------------------

  FairModule* cave= new MyCave("CAVE");
  cave->SetGeometryFileName("cave.geo");
  run->AddModule(cave);

  FairModule* magnet = new MyMagnet("Magnet");
  run->AddModule(magnet);

  FairModule* pipe = new MyPipe("Pipe");
  run->AddModule(pipe);
    
  FairDetector* NewDet = new NewDetector("TestDetector", kTRUE);
  run->AddModule(NewDet);

 // ------------------------------------------------------------------------


    // -----   Magnetic field   -------------------------------------------
    // Constant Field
    MyConstField  *fMagField = new MyConstField();
    fMagField->SetField(0., 20. ,0. ); // values are in kG
    fMagField->SetFieldRegion(-200, 200,-200, 200, -200, 200); // values are in cm
                          //  (xmin,xmax,ymin,ymax,zmin,zmax)
    run->SetField(fMagField);
    // --------------------------------------------------------------------

    
    
  // -----   Create PrimaryGenerator   --------------------------------------
  FairPrimaryGenerator* primGen = new FairPrimaryGenerator();
  
    // Add a box generator also to the run
    FairBoxGenerator* boxGen = new FairBoxGenerator(13, 5); // 13 = muon; 1 = multipl.
    boxGen->SetPRange(20,25); // GeV/c
    boxGen->SetPhiRange(0., 360.); // Azimuth angle range [degree]
    boxGen->SetThetaRange(0., 90.); // Polar angle in lab system range [degree]
    boxGen->SetXYZ(0., 0., 0.); // cm
    primGen->AddGenerator(boxGen);
  
    
    run->SetGenerator(primGen);
// ------------------------------------------------------------------------
 
  //---Store the visualiztion info of the tracks, this make the output file very large!!
  //--- Use it only to display but not for production!
  run->SetStoreTraj(kTRUE);

    
    
  // -----   Initialize simulation run   ------------------------------------
  run->Init();
  // ------------------------------------------------------------------------

  // -----   Runtime database   ---------------------------------------------

  Bool_t kParameterMerged = kTRUE;
  FairParRootFileIo* parOut = new FairParRootFileIo(kParameterMerged);
  parOut->open(parFile.Data());
  rtdb->setOutput(parOut);
  rtdb->saveOutput();
  rtdb->print();
  // ------------------------------------------------------------------------
   
  // -----   Start run   ----------------------------------------------------
   run->Run(nEvents);
    
  //You can export your ROOT geometry ot a separate file
  run->CreateGeometryFile("geofile_full.root");
  // ------------------------------------------------------------------------
  
  // -----   Finish   -------------------------------------------------------
  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
//.........这里部分代码省略.........

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;

}

//.........这里部分代码省略.........
      // ----- TRD reconstruction-----------------------------------------
      // Update of the values for the radiator F.U. 17.08.07
      Int_t trdNFoils = 130; // number of polyetylene foils
      Float_t trdDFoils = 0.0013; // thickness of 1 foil [cm]
      Float_t trdDGap = 0.02; // thickness of gap between foils [cm]
      Bool_t simpleTR = kTRUE; // use fast and simple version for TR production

      CbmTrdRadiator *radiator = new CbmTrdRadiator(simpleTR, trdNFoils, trdDFoils, trdDGap);

      if (trdHitProducerType == "smearing") {
         CbmTrdHitProducerSmearing* trdHitProd = new CbmTrdHitProducerSmearing(radiator);
         trdHitProd->SetUseDigiPar(false);
         run->AddTask(trdHitProd);
      } else if (trdHitProducerType == "digi") {
         CbmTrdDigitizer* trdDigitizer = new CbmTrdDigitizer(radiator);
         run->AddTask(trdDigitizer);

         CbmTrdHitProducerDigi* trdHitProd = new CbmTrdHitProducerDigi();
         run->AddTask(trdHitProd);
      } else if (trdHitProducerType == "clustering") {
         // ----- TRD clustering -----
         CbmTrdDigitizerPRF* trdClustering = new CbmTrdDigitizerPRF("TRD Clusterizer", "TRD task", radiator, false, true);
         run->AddTask(trdClustering);

         CbmTrdClusterFinderFast* trdClusterfindingfast = new CbmTrdClusterFinderFast(true, true, false, 5.0e-7);
         run->AddTask(trdClusterfindingfast);

         CbmTrdHitProducerCluster* trdClusterHitProducer = new CbmTrdHitProducerCluster();
         run->AddTask(trdClusterHitProducer);
         // ----- End TRD Clustering -----
      }
      // ------------------------------------------------------------------------
   }

   if (IsTof(parFile)) {
      // ------ TOF hits --------------------------------------------------------
      CbmTofHitProducerNew* tofHitProd = new CbmTofHitProducerNew("TOF HitProducerNew",iVerbose);
      tofHitProd->SetInitFromAscii(kFALSE);
      run->AddTask(tofHitProd);
      // ------------------------------------------------------------------------
   }

   // ------ Global track reconstruction -------------------------------------
   CbmLitFindGlobalTracks* finder = new CbmLitFindGlobalTracks();
   //CbmLitFindGlobalTracksParallel* finder = new CbmLitFindGlobalTracksParallel();
   // Tracking method to be used
   // "branch" - branching tracking
   // "nn" - nearest neighbor tracking
   // "nn_parallel" - nearest neighbor parallel tracking
   finder->SetTrackingType(std::string(globalTrackingType));

   // Hit-to-track merger method to be used
   // "nearest_hit" - assigns nearest hit to the track
   finder->SetMergerType("nearest_hit");

   run->AddTask(finder);

   if (IsTrd(parFile)) {
      CbmTrdMatchTracks* trdMatchTracks = new CbmTrdMatchTracks();
      run->AddTask(trdMatchTracks);
   }

   if (IsMuch(parFile)) {
      CbmMuchMatchTracks* muchMatchTracks = new CbmMuchMatchTracks();
      run->AddTask(muchMatchTracks);
   }

   // -----   Primary vertex finding   --------------------------------------
   CbmPrimaryVertexFinder* pvFinder = new CbmPVFinderKF();
   CbmFindPrimaryVertex* findVertex = new CbmFindPrimaryVertex(pvFinder);
   run->AddTask(findVertex);
   // -----------------------------------------------------------------------

	// -----  Parameter database   --------------------------------------------
	FairRuntimeDb* rtdb = run->GetRuntimeDb();
	FairParRootFileIo* parIo1 = new FairParRootFileIo();
	FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
	parIo1->open(parFile.Data());
	parIo2->open(parFileList, "in");
	rtdb->setFirstInput(parIo1);
	rtdb->setSecondInput(parIo2);
	rtdb->setOutput(parIo1);
	rtdb->saveOutput();
	// ------------------------------------------------------------------------

	// -----   Initialize and run   --------------------------------------------
	run->Init();
	run->Run(0, nEvents);
	// ------------------------------------------------------------------------

	// -----   Finish   -------------------------------------------------------
	timer.Stop();
 	cout << "Macro finished successfully." << endl;
   cout << "Test passed"<< endl;
   cout << " All ok " << endl;
	cout << "Output file is " << globalRecoFile << endl;
	cout << "Parameter file is " << parFile << endl;
	cout << "Real time " << timer.RealTime() << " s, CPU time "	<< timer.CpuTime() << " s" << endl;
	// ------------------------------------------------------------------------
}