void sim(Int_t nev=1) {

  gSystem->Load("liblhapdf");
  gSystem->Load("libEGPythia6");
  gSystem->Load("libpythia6");
  gSystem->Load("libAliPythia6");
  gSystem->Load("libhijing");
  gSystem->Load("libTHijing");
  gSystem->Load("libgeant321");

  if (gSystem->Getenv("EVENT"))
   nev = atoi(gSystem->Getenv("EVENT")) ;   
  
  AliSimulation simulator;
  simulator.SetMakeSDigits("TRD TOF PHOS HMPID EMCAL FMD ZDC PMD T0 VZERO");
  simulator.SetMakeDigitsFromHits("ITS TPC");
  simulator.SetWriteRawData("ALL","raw.root",kTRUE);

  simulator.SetDefaultStorage("local:///cvmfs/alice-ocdb.cern.ch/calibration/MC/Ideal");
  //simulator.SetDefaultStorage(Form("local://%s/OCDB", gSystem->pwd()));
  simulator.SetSpecificStorage("GRP/GRP/Data",
			       Form("local://%s",gSystem->pwd()));
  
  simulator.SetRunQA("ALL:ALL") ; 
  
  simulator.SetQARefDefaultStorage("local://$ALICE_ROOT/QAref") ;

  for (Int_t det = 0 ; det < AliQA::kNDET ; det++) {
    simulator.SetQACycles((AliQAv1::DETECTORINDEX_t)det, nev+1) ;
  }
  
  TStopwatch timer;
  timer.Start();
  simulator.Run(nev);
  timer.Stop();
  timer.Print();
}

void write(int n) {


  TRandom R;
  TStopwatch timer;




  TFile f1("mathcoreVectorIO_F.root","RECREATE");

  // create tree
  TTree t1("t1","Tree with new Float LorentzVector");

  XYZTVectorF *v1 = new XYZTVectorF();
  t1.Branch("LV branch","ROOT::Math::XYZTVectorF",&v1);

  timer.Start();
  for (int i = 0; i < n; ++i) {
     double Px = R.Gaus(0,10);
     double Py = R.Gaus(0,10);
     double Pz = R.Gaus(0,10);
     double E  = R.Gaus(100,10);
     //CylindricalEta4D<double> & c = v1->Coordinates();
     //c.SetValues(Px,pY,pZ,E);
     v1->SetCoordinates(Px,Py,Pz,E);
     t1.Fill();
  }

  f1.Write();
  timer.Stop();
  std::cout << " Time for new Float Vector " << timer.RealTime() << "  " << timer.CpuTime() << std::endl;

  t1.Print();

}

    void run_PPToGammaGammaFiles() {

    gROOT->LoadMacro("FlatTreeMaker_Delphes_PPToGammaGammaFiles_C.so");
  
    TChain* fChain = new TChain("Delphes");
    ifstream sourceFiles("PPToGammaGammaFiles.txt");
    char line[128];
    int  count = 0;
    cout<< "Adding files from PPToGammaGammaFiles to chain..."<< endl;
     while (sourceFiles >> line) {
        fChain->Add(line);
        ++count;
     }
    cout << count<<" files added!"<<endl;
    sourceFiles.close();
    TStopwatch timer;
    timer.Start();    
    fChain->Process("FlatTreeMaker_Delphes");

    cout << "\n\nDone!" << endl;
    cout << "CPU Time : " << timer.CpuTime() <<endl;
    cout << "RealTime : " << timer.RealTime() <<endl;                             
    cout <<"\n";
}

void makePlots() {

  gROOT->LoadMacro("analyze.C+");

  TStopwatch ts;
  ts.Start();

  TString input_ele = "ELE_FILE_TO_RUN";
  TString input_muon = "MUON_FILE_TO_RUN";
  bool addMC = true;
  int intLumi = 19712; // quote to 19.7

  double metCut = -1.;

  bool displayKStest = true;
  bool blinded = true;
  int nPhotons_req = 0;

  const int nChannels = 4;
  TString channels[nChannels] = {"ele_jjj", "ele_bjj",
				 "muon_jjj", "muon_bjj"};
  int nBtagReq[nChannels] = {0, 1,
			     0, 1};

  for(int i = 0; i < nChannels; i++) {
    if(i != 1 && i != 3) continue;
    if(i < 2) analyze(input_ele, addMC, i, intLumi, metCut, nPhotons_req, nBtagReq[i], displayKStest, blinded);
    else analyze(input_muon, addMC, i, intLumi, metCut, nPhotons_req, nBtagReq[i], displayKStest, blinded);
  }  

  ts.Stop();

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

}

void sim(Int_t nev=1) {
  AliSimulation simu;
  simu.SetMakeSDigits("TRD TOF PHOS HMPID  EMCAL MUON FMD PMD T0 ZDC VZERO");
  simu.SetMakeDigits ("TRD TOF PHOS HMPID  EMCAL MUON FMD PMD T0 ZDC VZERO");
  simu.SetMakeDigitsFromHits("ITS TPC");
  simu.SetWriteRawData("ALL","raw.root",kTRUE);
  simu.SetDefaultStorage("local://$ALICE_ROOT/OCDB");
  simu.SetSpecificStorage("GRP/GRP/Data",
			  Form("local://%s",gSystem->pwd()));

  simu.SetRunQA("ALL:ALL") ; 
  simu.SetQARefDefaultStorage("local://$ALICE_ROOT/OCDB") ;

  for (Int_t det = 0 ; det < AliQA::kNDET ; det++) {
    simu.SetQACycles(det, 2) ;
  }
  
  TStopwatch timer;
  timer.Start();
  simu.Run(nev);
  WriteXsection();
  timer.Stop();
  timer.Print();
}

int main(int argc, char **argv)
{
  TStopwatch reloj;
  reloj.Start();

  // split by ','
  string argStr = argv[1];
  vector<string> fileList;
  for (size_t i=0,n; i <= argStr.length(); i=n+1){
    n = argStr.find_first_of(',',i);
    if (n == string::npos) n = argStr.length();
    string tmp = argStr.substr(i,n-i);
    fileList.push_back(tmp);
  }

  Analysis1 o(fileList);
  o.EventsLoop();
  reloj.Stop();
  double tiempo = reloj.CpuTime();
  cout << "tiempo gastado en el calculo = " << tiempo << endl;

  return 0;

}

void sim(Int_t nev=1) {

  gSystem->Exec(" rm itsSegmentations.root ");
  AliSimulation simulator;
  //  simulator.SetMakeSDigits("");

  //  simulator.SetMakeDigits("");

  simulator.SetDefaultStorage("local://$ALICE_ROOT/OCDB");
  simulator.SetSpecificStorage("GRP/GRP/Data",
			       Form("local://%s",gSystem->pwd()));
  simulator.SetSpecificStorage("ITS/Align/Data",
			       Form("local://%s",gSystem->pwd()));
  simulator.SetSpecificStorage("ITS/Calib/SimuParam",
			       Form("local://%s",gSystem->pwd()));
  simulator.SetRunHLT("");
  simulator.SetRunQA(":");

  TStopwatch timer;
  timer.Start();
  simulator.Run(nev);
  timer.Stop();
  timer.Print();
}

void  testIntegPerf(double x1, double x2, int n = 100000){


   std::cout << "\n\n***************************************************************\n";
   std::cout << "Test integration performances in interval [ " << x1 << " , " << x2 << " ]\n\n";

  TStopwatch timer;

  double dx = (x2-x1)/double(n);

  //ROOT::Math::Functor1D<ROOT::Math::IGenFunction> f1(& TMath::BreitWigner);
  ROOT::Math::WrappedFunction<> f1(func);

  timer.Start();
  ROOT::Math::Integrator ig(f1 );
  double s1 = 0.0;
  nc = 0;
  for (int i = 0; i < n; ++i) {
     double x = x1 + dx*i;
     s1+= ig.Integral(x1,x);
  }
  timer.Stop();
  std::cout << "Time using ROOT::Math::Integrator        :\t" << timer.RealTime() << std::endl;
  std::cout << "Number of function calls = " << nc/n << std::endl;
  int pr = std::cout.precision(18);  std::cout << s1 << std::endl;  std::cout.precision(pr);



  //TF1 *fBW = new TF1("fBW","TMath::BreitWigner(x)",x1, x2);  //  this is faster but cannot measure number of function calls
  TF1 *fBW = new TF1("fBW",func2,x1, x2,0);

  timer.Start();
  nc = 0;
  double s2 = 0;
  for (int i = 0; i < n; ++i) {
     double x = x1 + dx*i;
     s2+= fBW->Integral(x1,x );
  }
  timer.Stop();
  std::cout << "Time using TF1::Integral :\t\t\t" << timer.RealTime() << std::endl;
  std::cout << "Number of function calls = " << nc/n << std::endl;
  pr = std::cout.precision(18);  std::cout << s1 << std::endl;  std::cout.precision(pr);


}

// test using UNURAN string interface
void testStringAPI() {

   TH1D * h1 = new TH1D("h1G","gaussian distribution from Unuran",100,-10,10);
   TH1D * h2 = new TH1D("h2G","gaussian distribution from TRandom",100,-10,10);

   cout << "\nTest using UNURAN string API \n\n";


   TUnuran unr; 
   if (! unr.Init( "normal()", "method=arou") ) {
      cout << "Error initializing unuran" << endl;
      return;
   }

   int n = NGEN;
   TStopwatch w; 
   w.Start(); 

   for (int i = 0; i < n; ++i) {
       double x = unr.Sample(); 
       h1->Fill(  x ); 
   }

   w.Stop(); 
   cout << "Time using Unuran method " << unr.MethodName() << "\t=\t " << w.CpuTime() << endl;


   // use TRandom::Gaus
   w.Start();
   for (int i = 0; i < n; ++i) {
      double x = gRandom->Gaus(0,1); 
       h2->Fill(  x ); 
   }

   w.Stop(); 
   cout << "Time using TRandom::Gaus  \t=\t " << w.CpuTime() << endl;

   assert(c1 != 0);
   c1->cd(++izone);
   h1->Draw();
   c1->cd(++izone);
   h2->Draw();

}

void testDiscDistr() { 

   cout << "\nTest Discrete distributions\n\n";

   TH1D * h1 = new TH1D("h1PS","Unuran Poisson prob",20,0,20);
   TH1D * h2 = new TH1D("h2PS","Poisson dist from TRandom",20,0,20);

   double mu = 5; 

   TF1 * f = new TF1("fps",poisson,1,0,1);
   f->SetParameter(0,mu);

   TUnuranDiscrDist dist2 = TUnuranDiscrDist(f);
   TUnuran unr;

   // dari method (needs also the mode and pmf sum)
   dist2.SetMode(int(mu) );
   dist2.SetProbSum(1.0);
   bool ret = unr.Init(dist2,"dari");
   if (!ret) return;

   TStopwatch w; 
   w.Start(); 

   int n = NGEN;
   for (int i = 0; i < n; ++i) {
      int k = unr.SampleDiscr(); 
      h1->Fill( double(k) ); 
   }

   w.Stop(); 
   cout << "Time using Unuran method " << unr.MethodName() << "\t=\t\t " << w.CpuTime() << endl;

   w.Start();
   for (int i = 0; i < n; ++i) {
      h2->Fill(  gRandom->Poisson(mu) ); 
   }
   cout << "Time using TRandom::Poisson " << "\t=\t\t " << w.CpuTime() << endl;

   c1->cd(++izone);
   h1->SetMarkerStyle(20);
   h1->Draw("E");
   h2->Draw("same");
   
   std::cout << " chi2 test of UNURAN vs TRandom generated histograms:  " << std::endl;
   h1->Chi2Test(h2,"UUP");

}

void generate( R & r, TH1D * h) { 

  TStopwatch w; 

  r.SetSeed(0);
  //r.SetSeed(int(std::pow(2.0,28)));
  int m = NLOOP;
  int n = NEVT;
  for (int j = 0; j < m; ++j) { 

    //std::cout << r.GetSeed() << "   "; 

    w.Start();
//     if ( n < 40000000) iseed = std::rand();
//     iseed = 0;
    //TRandom3 r3(0);
    //r.SetSeed( 0 ); // generate random seeds
    //TRandom3 r3(0); 
    //r.SetSeed (static_cast<UInt_t> (4294967296.*r3.Rndm()) );

  // estimate PI
    double n1=0; 
    double rn[2000];
    double x; 
    double y; 
    for (int ievt = 0; ievt < n; ievt+=1000 ) { 
      r.RndmArray(2000,rn);
      for (int i=0; i < 1000; i++) { 
	x=rn[2*i];
	y=rn[2*i+1];
	if ( ( x*x + y*y ) <= 1.0 ) n1++;
      }
    }
    double piEstimate = 4.0 * double(n1)/double(n);
    double delta = piEstimate-PI; 
    h->Fill(delta); 
  }

  w.Stop();
  std::cout << std::endl; 
  std::cout << "Random:  " << typeid(r).name() 
	    << "\n\tTime = " << w.RealTime() << "  " << w.CpuTime() << std::endl;   
  std::cout << "Time/call:  " << w.CpuTime()/(2*n)*1.0E9 << std::endl; 
}

//.........这里部分代码省略.........
         }
         else { 
            Error("StandardHypoTestInvDemo","Model %s has no valid poi",modelBName);
            return 0;
         }         
      }
   }

   // check model  has global observables when there are nuisance pdf
   // for the hybrid case the globobs are not needed
   if (type != 1 ) { 
      bool hasNuisParam = (sbModel->GetNuisanceParameters() && sbModel->GetNuisanceParameters()->getSize() > 0);
      bool hasGlobalObs = (sbModel->GetGlobalObservables() && sbModel->GetGlobalObservables()->getSize() > 0);
      if (hasNuisParam && !hasGlobalObs ) {  
         // try to see if model has nuisance parameters first 
         RooAbsPdf * constrPdf = RooStats::MakeNuisancePdf(*sbModel,"nuisanceConstraintPdf_sbmodel");
         if (constrPdf) { 
            Warning("StandardHypoTestInvDemo","Model %s has nuisance parameters but no global observables associated",sbModel->GetName());
            Warning("StandardHypoTestInvDemo","\tThe effect of the nuisance parameters will not be treated correctly ");
         }
      }
   }


  
   // 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.....");

void inflateTree(const char *name = "h42",
                 const char *in = "root://eospps.cern.ch///eos/ppsscratch/test/h1big.root",
                 const char *out = "/tmp/h1big.root",  Int_t fact = 1)
{
	TStopwatch sw;
	sw.Start();
	
   // Get the input tree from the input file
   TFile *fin = TFile::Open(in);
   if (!fin || fin->IsZombie()) {
      Printf("inflateTree", "could not open input file: %s", in);
      return;
   }
   TTree *tin = (TTree *) fin->Get(name);
   if (!tin) {
      Printf("inflateTree", "could not find tree '%s' in %s", name, in);
      delete fin;
      return;
   }
   Long64_t nin = tin->GetEntriesFast();
   Printf("Input tree '%s' has %lld entries", name, nin);
   // Create output file
   TFile *fout = TFile::Open(out, "RECREATE", 0, 1);
   if (!fout || fout->IsZombie()) {
      Printf("inflateTree", "could not open input file: %s", in);
      delete fin;
      return;
   }
   // Clone the header of the initial tree
   TTree *tout= (TTree *) tin->CloneTree(0);
   tout->SetMaxTreeSize(19000000000);


   // Duplicate all entries once
#if 0
   Int_t nc = fact;
   while (nc--) {
      Printf("Writing copy %d ...", fact - nc);
      tout->CopyEntries(tin);
   }
#else
   for (Long64_t i = 0; i < nin; ++i) {
      if (tin->LoadTree(i) < 0) {
         break;
      }
      tin->GetEntry(i);
      Int_t nc = fact;
      while (nc--) {
         tout->Fill();
      }
      if (i > 0 && !(i%1000)) {
         Printf("%d copies of %lld entries filled ...", fact, i);
      }
   }
#endif
   // Finalize the writing out
   tout->Write();
   
   // print perf stats
   

    sw.Stop();
    std::cout << "Drawing. Realtime: " <<      sw.RealTime()  << std::endl;
    std::cout << "Drawing. Cputime : " <<      sw.CpuTime()  << std::endl;
	tin->PrintCacheStats();   
   
   
   // Close the files
   fout->Close();
   fin->Close();
   // Cleanup
   delete fout;
   delete fin;
}

void r3blandreco(Int_t nNeutrons, Int_t beamE, Int_t Erel)
{
  Int_t d;
  if(Erel == 100){
    d = 35;
  }
  else{
    d = 14;
  }
  
  // ----- Files ---------------------------------------------------------------
  char strDir[] = ".";
  char str[100];
  char str2[100];
  sprintf(str, "%1dAMeV.%1dn.%1dkeV.%1dm.root", beamE,nNeutrons, Erel, d);
  sprintf(str2, "%1dAMeV.%1dkeV.%1dm", beamE, Erel, d);
  TString inFile  = TString(strDir) + "/r3bsim." + TString(str);
  TString digiFile  = TString(strDir) + "/r3bcalibr." + TString(str);
  TString parFile  = TString(strDir) + "/r3bpar." + TString(str);
  TString calibrFile  = TString(strDir) + "/r3bcalibr." + TString(str2) + ".txt";
  TString outFile  = TString(strDir) + "/r3breco." + TString(str);
  // ---------------------------------------------------------------------------



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



  // ----- Digitization --------------------------------------------------------
  FairRunAna *fRun= new FairRunAna();
  fRun->SetInputFile(inFile);
  fRun->AddFriend(digiFile);
  fRun->SetOutputFile(outFile);
  // ---------------------------------------------------------------------------


  // ---------------------------------------------------------------------------
  Double_t beamEnergy;
  Double_t beamBeta;
  if(200 == beamE) {
    beamEnergy=200.;
    beamBeta=0.5676881;
  } else if(600 == beamE) {
    beamEnergy=600.;
    beamBeta=0.7937626;
  } else if(1000 == beamE) {
    beamEnergy=1000.;
    beamBeta=0.8760237;
  }
  // ---------------------------------------------------------------------------
  
  
  
  // ----- Connect the Tracking Task -------------------------------------------
  R3BNeutronTracker2D* tracker  = new R3BNeutronTracker2D();
  tracker->UseBeam(beamEnergy, beamBeta);
  tracker->ReadCalibrFile(calibrFile.Data());
  fRun->AddTask(tracker);
  // ---------------------------------------------------------------------------



  // ----- Runtime DataBase info -----------------------------------------------
  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
  FairParRootFileIo*  parIo1 = new FairParRootFileIo();
  parIo1->open(parFile.Data());
  rtdb->setFirstInput(parIo1);
  rtdb->setOutput(parIo1);
  rtdb->saveOutput();
  // ---------------------------------------------------------------------------



  // ----- Number of events to process -----------------------------------------
  Int_t nEvents = 10000;
  // ---------------------------------------------------------------------------


  
  // ----- Intialise and run ---------------------------------------------------
  fRun->Init();
  fRun->Run(0, 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 writen:  "    << outFile << endl;
  cout << "Parameter file writen " << parFile << endl;
  cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
  cout << endl;
//.........这里部分代码省略.........

void run(const Char_t *files=NULL, Bool_t mc=kFALSE, Bool_t tpid=kTRUE,  Bool_t tchg=kFALSE,  Bool_t tpp=kTRUE, Long64_t nev=1234567890, Long64_t first = 0)
{
  TStopwatch timer;
  timer.Start();

  // VERY GENERAL SETTINGS
  //AliLog::SetGlobalLogLevel(AliLog::kError);
  if(gSystem->Load("libANALYSIS.so")<0) return;
  if(gSystem->Load("libANALYSISalice.so")<0) return;
  if(gSystem->Load("libTender.so")<0) return;
  if(gSystem->Load("libTenderSupplies.so")<0) return;
//   if(gSystem->Load("libMES.so")<0) return;
    if(gSystem->Load("libPWGLFspectra.so")<0) return;

  // DEFINE DATA CHAIN
  TChain *chain = NULL;
  if(!files) chain = MakeChainLST();
  else chain = MakeChainLST(files);

  if(!chain) return;
  chain->Lookup();
  chain->GetListOfFiles()->Print();
  Long64_t nfound=chain->GetEntries();
  printf("\tENTRIES FOUND [%lli] REQUESTED [%lli]\n", nfound, nev>nfound?nfound:nev);

  // BUILD ANALYSIS MANAGER
  AliAnalysisManager *mgr = new AliAnalysisManager("Multiplicity and Event Shape");
  AliESDInputHandler *esdH = new AliESDInputHandler();
  AliMCEventHandler *mcH(NULL);
  mgr->SetInputEventHandler(esdH);
  if(mc) mgr->SetMCtruthEventHandler(mcH = new AliMCEventHandler());
  //mgr->SetDebugLevel(10);
  mgr->SetSkipTerminate(kTRUE);

  // LOAD tasks
  // *******************  PID response  ******************
  gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPIDResponse.C");
  if(!mc) AddTaskPIDResponse();
  else AddTaskPIDResponse(kTRUE,kTRUE,kTRUE,2);

  // *******************  Tenders  ***********************
  AliTender *aliTender(NULL);
  gROOT->LoadMacro("$ALICE_PHYSICS/TENDER/TenderSupplies/AddTaskTender.C");
  if(!mc){    // for DATA
    aliTender = (AliTender*)AddTaskTender(!mc, kTRUE, kTRUE, kTRUE, kTRUE, kFALSE, kTRUE, kFALSE, kFALSE);
       // (useV0, useTPC,  !!! useTOF=kFALSE for MC !!!, useTRD, usePID, useVTX, useT0, useEmc, usePtFix)
  } else {  // for MC
    aliTender = (AliTender*)AddTaskTender(!mc, kTRUE, kFALSE, kTRUE, kTRUE, kTRUE, kTRUE, kFALSE, kFALSE);  // (useV0, useTPC,  !!! useTOF=kFALSE for MC !!!, useTRD, usePID, useVTX, useT0, useEmc, usePtFix)
  }
  aliTender->SetHandleOCDB(kTRUE);
  //aliTender->SetDefaultCDBStorage(Form("alien://folder=/alice/data/2010/OCDB?cacheFolder=%s/local", gSystem->ExpandPathName("$HOME")));
  // aliTender->SetDefaultCDBStorage(Form("local://%s/local/alice/data/2010/OCDB", gSystem->ExpandPathName("$HOME")));

// *******************  Physics Selection  *************
  gROOT->LoadMacro("$ALICE_PHYSICS/OADB/macros/AddTaskPhysicsSelection.C");
  AliPhysicsSelectionTask *physSelTask = AddTaskPhysicsSelection(mc); // 0 = real data; 1 = MC

  // *******************  MES Tender  ******************
  gROOT->LoadMacro("$ALICE_PHYSICS/PWGLF/SPECTRA/MultEvShape/AddMEStender.C");
  AddMEStender(mc);

  // *******************  MES PID task  ******************
  if(tpid){
    gROOT->LoadMacro("$ALICE_PHYSICS/PWGLF/SPECTRA/MultEvShape/AddMESpidTask.C");
	AddMESpidTask(mc);
  }
//

//   // *******************  MES CHG task  ******************
  if(tchg){
    gROOT->LoadMacro("$ALICE_PHYSICS/PWGLF/SPECTRA/MultEvShape/AddMESchgTask.C");
    AddMESchgTask(mc);
  }
//
//   // *******************  MES ppCol task  ******************
  if(tpp){
    gROOT->LoadMacro("$ALICE_PHYSICS/PWGLF/SPECTRA/MultEvShape/AddMESppColTask.C");
    AddMESppColTask(mc);
  }


  if (!mgr->InitAnalysis()) return;
  mgr->PrintStatus();
  mgr->StartAnalysis("local", chain, nev, first);
  timer.Stop();
  timer.Print();
  // verbosity
  printf("\tCLEANING TASK LIST:\n");
  mgr->GetTasks()->Delete();

  if(mcH) delete mcH;
  delete esdH;
  delete chain;
}

//.........这里部分代码省略.........
//std::vector<double> *myleptoncharge=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(Bmodification((*mycvsdiscriminant)[0],(*myjetmatchinginfo)[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( Bmodification((*mycvsdiscriminant)[0],(*myjetmatchinginfo)[0]),finalweight );
if (samples_[idx]=="WPHJET.root")outsidewphjet=outsidewphjet+finalweight;}


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

int quasirandom(int n = 10000, int skip = 0) {

   TH2D * h0 = new TH2D("h0","Pseudo-random Sequence",200,0,1,200,0,1);
   TH2D * h1 = new TH2D("h1","Sobol Sequence",200,0,1,200,0,1);
   TH2D * h2 = new TH2D("h2","Niederrer Sequence",200,0,1,200,0,1);

   RandomMT         r0;
   // quasi random numbers need to be created giving the dimension of the sequence
   // in this case we generate 2-d sequence

   QuasiRandomSobol r1(2);
   QuasiRandomNiederreiter r2(2);

   // generate n random points

   double x[2];
   TStopwatch w; w.Start();
   for (int i = 0; i < n; ++i)  {
      r0.RndmArray(2,x);
      h0->Fill(x[0],x[1]);
   }
   std::cout << "Time for gRandom ";
   w.Print();

   w.Start();
   if( skip>0) r1.Skip(skip);
   for (int i = 0; i < n; ++i)  {
      r1.Next(x);
      h1->Fill(x[0],x[1]);
   }
   std::cout << "Time for Sobol ";
   w.Print();

   w.Start();
   if( skip>0) r2.Skip(skip);
   for (int i = 0; i < n; ++i)  {
      r2.Next(x);
      h2->Fill(x[0],x[1]);
   }
   std::cout << "Time for Niederreiter ";
   w.Print();

   TCanvas * c1 = new TCanvas("c1","Random sequence",600,1200);
   c1->Divide(1,3);
   c1->cd(1);
   h0->Draw("COLZ");
   c1->cd(2);

   // check uniformity
   h1->Draw("COLZ");
   c1->cd(3);
   h2->Draw("COLZ");
   gPad->Update();

   // test number of empty bins

   int nzerobins0 = 0;
   int nzerobins1 = 0;
   int nzerobins2 = 0;
   for (int i = 1; i <= h1->GetNbinsX(); ++i) {
      for (int j = 1; j <= h1->GetNbinsY(); ++j) {
         if (h0->GetBinContent(i,j) == 0 ) nzerobins0++;
         if (h1->GetBinContent(i,j) == 0 ) nzerobins1++;
         if (h2->GetBinContent(i,j) == 0 ) nzerobins2++;
      }
   }

   std::cout << "number of empty bins for pseudo-random = " << nzerobins0 << std::endl;
   std::cout << "number of empty bins for " << r1.Name() << "\t= " << nzerobins1 << std::endl;
   std::cout << "number of empty bins for " << r2.Name() << "\t= " << nzerobins2 << std::endl;

   int iret = 0;
   if (nzerobins1 >= nzerobins0 ) iret += 1;
   if (nzerobins2 >= nzerobins0 ) iret += 2;
   return iret;

}

RooGaussian fitZToMuMuGammaMassUnbinned(const char *filename = "ZToMuMuGammaMass.txt",
  const char* plotOpt = "NEU",
  const int nbins = 25)
{
  gROOT->ProcessLine(".L tdrstyle.C");
  setTDRStyle();
  gStyle->SetPadRightMargin(0.05);

  double minMass = 60;
  double maxMass = 120;
  RooRealVar  mass("mass","M(#mu#mu#gamma})", minMass, maxMass,"GeV/c^{2}");

  // Read data set

  RooDataSet *data = RooDataSet::read(filename,RooArgSet(mass));
//   RooDataSet *dataB = RooDataSet::read(filenameB,RooArgSet(mass));

// Build p.d.f.

////////////////////////////////////////////////
//             Parameters                     //
////////////////////////////////////////////////

//  Signal p.d.f. parameters
//  Parameters for a Gaussian and a Crystal Ball Lineshape
  RooRealVar  m0   ("m_{0}", "Bias", 91.19, minMass, maxMass,"GeV/c^{2}");
  RooRealVar  sigma("#sigma","Width", 3.0,1.0,10.0,"GeV/c^{2}");
  RooRealVar  cut  ("#alpha","Cut", 0.6,0.6,2.0);
  RooRealVar  power("power","Power", 10.0, 0.5, 20.0);


//  Background p.d.f. parameters
//  Parameters for a polynomial lineshape
  RooRealVar  c0("c_{0}", "c0", 0., -10, 10);
  RooRealVar  c1("c_{1}", "c1", 0., -100, 0);
  RooRealVar  c2("c_{2}", "c2", 0., -100, 100);
//  c0.setConstant();

// fraction of signal
//  RooRealVar  frac("frac", "Signal Fraction", 0.1,0.,0.3.);
  RooRealVar  nsig("N_{S}", "#signal events", 9000, 0.,10000.);
  RooRealVar  nbkg("N_{B}", "#background events", 1000,2,10000.);



////////////////////////////////////////////////
//               P.D.F.s                      //
////////////////////////////////////////////////

// Di-photon mass signal p.d.f.
  RooGaussian    signal("signal", "A  Gaussian Lineshape", mass, m0, sigma);
  // RooCBShape     signal("signal", "A  Crystal Ball Lineshape", mass, m0,sigma, cut, power);

// Di-photon mass background  p.d.f.
  RooPolynomial bg("bg", "Backgroung Distribution", mass, RooArgList(c0,c1));

// Di-photon mass model p.d.f.
//  RooAddPdf      model("model", "Di-photon mass model", signal, bg, frac);
  RooAddPdf      model("model", "Di-photon mass model", RooArgList(signal, bg), RooArgList(nsig, nbkg));


  TStopwatch t ;
  t.Start() ;
//   model->fitTo(*data,FitOptions("mh"),Optimize(0),Timer(1));
  signal->fitTo(*data,FitOptions("mh"),Optimize(0),Timer(1));

  t.Print() ;

  c = new TCanvas("c","Unbinned Invariant Mass Fit", 0,0,800,600);
// Plot the fit results
  RooPlot* plot = mass.frame(Range(minMass,maxMass),Bins(nbins));

// Plot 1
//   dataB->plotOn(plot, MarkerColor(kRed), LineColor(kRed));
  data->plotOn(plot);
//   model.plotOn(plot);
  model.plotOn(plot);
  //model.paramOn(plot, Format(plotOpt, AutoPrecision(1)), Parameters(RooArgSet(nsig, nbkg, m0, sigma)));
  model.paramOn(plot, Format(plotOpt, AutoPrecision(1)), Parameters(RooArgSet(m0, sigma)));

  /// model.plotOn(plot, Components("signal"), LineStyle(kDashed), LineColor(kRed));

//   model.plotOn(plot, Components("bg"), LineStyle(kDashed), LineColor(kRed));


  plot->Draw();

  TLatex *   tex = new TLatex(0.2,0.8,"CMS preliminary");
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetLineWidth(2);
  tex->Draw();
  tex->DrawLatex(0.2, 0.725, "7 TeV Data, L = 258 pb^{-1}");

  float fsig_peak = NormalizedIntegral(signal,
                      mass,
                      m0.getVal() - 2.5*sigma.getVal(),
                      m0.getVal() + 2.5*sigma.getVal()
                    );

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

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\&quo