C++ TH1::Fill方法代码示例

void parallelMergeClient() 
{
   // Client program which creates and fills 2 histograms and a TTree. 
   // Every 1000000 fills the histograms and TTree is send to the server which displays the histogram.
   //
   // To run this demo do the following:
   //   - Open at least 2 windows
   //   - Start ROOT in the first windows
   //   - Execute in the first window: .x fastMergeServer.C
   //   - Execute in the other windows: root.exe -b -l -q .x treeClient.C
   //     (You can put it in the background if wanted).
   // If you want to run the hserv.C on a different host, just change
   // "localhost" in the TSocket ctor below to the desired hostname.
   //
   //Author: Fons Rademakers, Philippe Canal
   
   gBenchmark->Start("treeClient");

   TParallelMergingFile *file = (TParallelMergingFile*)TFile::Open("mergedClient.root?pmerge=localhost:1095","RECREATE");
   
   file->Write();
   file->UploadAndReset();       // We do this early to get assigned an index.
   UInt_t idx = file->fServerIdx; // This works on in ACLiC.

   TH1 *hpx;
   if (idx == 0) {
      // Create the histogram
      hpx = new TH1F("hpx","This is the px distribution",100,-4,4);
      hpx->SetFillColor(48);  // set nice fillcolor
   } else {
      hpx = new TH2F("hpxpy","py vs px",40,-4,4,40,-4,4);
   }
   Float_t px, py;
   TTree *tree = new TTree("tree","tree");
   tree->SetAutoFlush(4000000);
   tree->Branch("px",&px);
   tree->Branch("py",&py);
 
   // Fill histogram randomly
   gRandom->SetSeed();
   const int kUPDATE = 1000000;
   for (int i = 0; i < 25000000; ) {
      gRandom->Rannor(px,py);
      if (idx%2 == 0)
         hpx->Fill(px);
      else
         hpx->Fill(px,py);
      tree->Fill();
      ++i;
      if (i && (i%kUPDATE) == 0) {
         file->Write();
      }
   }
   file->Write();
   delete file;

   gBenchmark->Show("treeClient");
}

// PValue for finding a local p-value as observed in 'bin' or worse 
void ToyExperiments::printGlobalPValueOfLocalFluctuation(unsigned int bin, unsigned int nExperiments) const {
  std::cout << "Determining global p-value for observed fluctuation in bin " << bin << " from " << nExperiments << " toy experiments ...  " << std::flush;

  // Find the predicted yields that correspond
  // to the local p-value 'localPValue'
  std::vector<unsigned int> limitYields = yields(localPValue(bin,observedYields_.at(bin)));

  TH1* hIsAbovePValue = new TH1D("hIsAbovePValue","",2,0,2);

  const double minCorr = findMinValidRandomNumberForCorrelatedUncertainties();

  for(unsigned int p = 0; p < nExperiments; ++p) {
    bool isAbovePValue = false;
    double rCorr = rand_->Gaus(0.,1.);
    while( rCorr <= minCorr ) {
      rCorr = rand_->Gaus(0.,1.);
    }
    for(unsigned int b = 0; b < Parameters::nBins(); ++b) {
      double prediction = -1.;
      bool negativePrediction = true;
      while( negativePrediction ) {
	double rUncorr = rand_->Gaus(0.,1.);
	double uncorrVar = rUncorr * uncorrelatedUncerts_.at(b);
	double corrVar   = rCorr   * correlatedUncerts_.at(b);
	prediction = meanPredictions_.at(b) + uncorrVar + corrVar;
	if( prediction >= 0. ) {
	  negativePrediction = false;
	}
      }
      double predictedYield = rand_->Poisson(prediction);
      if( predictedYield >= limitYields.at(b) ) {
	isAbovePValue = true;
	break;
      }      
    }
    if( isAbovePValue ) {
      hIsAbovePValue->Fill(1);
    } else {
      hIsAbovePValue->Fill(0);
    }
  }
  std::cout << "ok" << std::endl;

  double lpv      = localPValue(bin,observedYields_.at(bin));
  double gpUncorr = 1. - pow(1.-lpv,Parameters::nBins());
  double gpCorr   = hIsAbovePValue->Integral(2,2)/hIsAbovePValue->Integral(1,2);

  std::cout << "\n\n----- Global p-value for observed fluctuation in bin " << bin << " -----" << std::endl;
  std::cout << "  local p-value                           : " << lpv << " (" << TMath::NormQuantile(1.-lpv) << "sig)" << std::endl;
  std::cout << "  global p-value (without correlations)   : " << gpUncorr  << " (" << TMath::NormQuantile(1.-gpUncorr) << "sig)" << std::endl;
  std::cout << "  global p-value (including correlations) : " << gpCorr  << " (" << TMath::NormQuantile(1.-gpCorr) << "sig)" << std::endl;
}

TH1* fillHistogram(TTree* testTree, const std::string& varName, const std::string& selection, const std::string& frWeightName,
		   const std::string& histogramName, unsigned numBinsX, double xMin, double xMax)
{
  std::cout << "<fillHistogram>:" << std::endl;
  std::cout << " testTree = " << testTree << std::endl;
  std::cout << " varName = " << varName << std::endl;
  std::cout << " selection = " << selection << std::endl;
  std::cout << " frWeightName = " << frWeightName << std::endl;
  std::cout << " histogramName = " << histogramName << std::endl;

  TH1* histogram = new TH1F(histogramName.data(), histogramName.data(), numBinsX, xMin, xMax);
  histogram->Sumw2();

  TFile* dummyOutputFile = new TFile("dummyOutputFile.root", "RECREATE");

  TTree* selTree = ( selection != "" ) ? testTree->CopyTree(selection.data()) : testTree;
  std::cout << " selTree = " << selTree << std::endl;

  Float_t eventWeight = 1.;
  selTree->SetBranchAddress("weight", &eventWeight);

  Float_t var = 0.;
  selTree->SetBranchAddress(varName.data(), &var);
  
  Float_t frWeight = 1.;
  if ( frWeightName != "" ) {
    std::cout << "--> setting branch-address of fake-rate weight..." << std::endl;
    selTree->SetBranchAddress(frWeightName.data(), &frWeight);
  }

  int numEntries = selTree->GetEntries();
  std::cout << "--> numEntries = " << numEntries << std::endl;
  //if ( maxEntriesToProcess != -1 ) numEntries = TMath::Min(numEntries, maxEntriesToProcess);
  for ( int iEntry = 0 ; iEntry < numEntries; ++iEntry ) {
    selTree->GetEvent(iEntry);  

    //std::cout << "iEntry = " << iEntry << ": var = " << var << ", frWeight = " << frWeight << std::endl;

    if ( frWeightName != "" ) {
      if ( TMath::Abs(frWeight) < 1. ) // some entries have weight O(-100)
                                       // --> indication of technical problem with k-NearestNeighbour tree ?
	histogram->Fill(var, frWeight*eventWeight);
    } else {
      histogram->Fill(var, eventWeight);
    }
  }

  delete dummyOutputFile;

  return histogram;
}

int main (int argc, char *argv[])
{ 
	LineStream ls;
	char *line;
	char *pos;
	Stringa buffer;

	if (argc != 2) {
		usage ("%s <file.intraOffsets>");
	}

	TH1 *his = new TH1D ("","Intra-read distribution",1000,0,1000);
	TCanvas *canv = new TCanvas("","canvas",1200,400);
	ls = ls_createFromFile (argv[1]);
	while (line = ls_nextLine (ls)) {
		his->Fill (atoi (line));
	}
	ls_destroy (ls);
	his->Draw();
	his->GetXaxis()->SetLabelSize (0.04);
	his->GetYaxis()->SetLabelSize (0.04);
	buffer = stringCreate (100);
	pos = strchr (argv[1],'.');
	if (pos == NULL) {
		die ("Expected <file.intraOffsets>: %s",argv[1]);
	}
	*pos = '\0';
	stringPrintf (buffer,"%s_intraDistribution.jpg",argv[1]);
	canv->Print (string (buffer),"jpg");
	stringDestroy (buffer);
	return 0;
}

void Kin2Txt() {
  AliRunLoader* rl = AliRunLoader::Open("galice.root");
  rl->LoadKinematics();
  rl->LoadHeader();

  TH1* hM  = new TH1D("hM",  "TreeK;M#(){#pi^{+}#pi^{-}} #(){GeV/#it{c}^{2}}", 100, 0., 2.);
  TH1* hPt = new TH1D("hPt", "TreeK;P_{T}#(){#pi^{+}#pi^{-}} #(){GeV/#it{c}}", 100, 0., 1.);
  TH1* hY  = new TH1D("hY",  "TreeK;Y#(){#pi^{+}#pi^{-}}",                     160,-8., 8.);
  std::ofstream ofs("rho0.txt");
  
  AliStack *stack = NULL;
  TParticle *part = NULL;
  TLorentzVector v[2], vSum;
  for (Int_t i=0, n(rl->GetNumberOfEvents()); i<n; ++i) {
    rl->GetEvent(i);
    stack = rl->Stack();

    Int_t nPrimary(0);
    for (Int_t j(0), m(stack->GetNtrack()); j<m; ++j) {
      part = stack->Particle(j);
      if (part->IsPrimary())
	part->Momentum(v[nPrimary++]);
    }
    if (nPrimary != 2) {
      Printf("Error: nPrimary=%d != 2", nPrimary);
      continue;
    }
    vSum = v[0] + v[1];
    hM->Fill(vSum.M());
    hPt->Fill(vSum.Perp());
    hY->Fill(vSum.Rapidity());
    ofs << std::fixed << std::setprecision(4)
	<< vSum.M() << " " << vSum.Perp() << " " << vSum.Rapidity() << " "
	<< v[0].Eta() << " " << v[0].Px() << " " << v[0].Py() << " " << v[0].Pz() << " "
	<< v[1].Eta() << " " << v[1].Px() << " " << v[1].Py() << " " << v[1].Pz()
	<< std::endl;
  }
  hM->Draw();
  c1->SaveAs("TreeK.pdf(");
  hPt->Draw();
  c1->SaveAs("TreeK.pdf");
  hY->Draw();
  c1->SaveAs("TreeK.pdf)");
}

void analyzer()
{

  TString processName = "ZJets";

  TFile* f = TFile::Open(Form("hist_%s.root", processName.Data()), "recreate");

  // Create chain of root trees
  TChain chain("DelphesMA5tune");
  //kisti
  //chain.Add("/cms/home/tjkim/fcnc/sample/ZToLL50-0Jet_sm-no_masses/events_*.root");
  //hep
  chain.Add("/Users/tjkim/Work/fcnc/samples/ZToLL50-0Jet_sm-no_masses/events_*.root");  

  // Create object of class ExRootTreeReader
  ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
  Long64_t numberOfEntries = treeReader->GetEntries();
  
  // Get pointers to branches used in this analysis
  TClonesArray *branchMuon = treeReader->UseBranch("DelphesMA5tuneMuon");
  
  // Book histograms
  TH1 *histDiMuonMass = new TH1F("dimuon_mass","Di-Muon Invariant Mass (GeV)",100, 50, 150);
 
  // Loop over all events

  for(Int_t entry = 0; entry < numberOfEntries; ++entry)
  {
    // Load selected branches with data from specified event
    treeReader->ReadEntry(entry);
 
    int nmuon = 0; 
    for( int i = 0; i < branchMuon->GetEntries(); i++)
    {
      Muon *muon = (Muon*) branchMuon->At(i);
      if( muon->PT <= 20 || abs(muon->Eta) >= 2.4 ) continue;
      nmuon = nmuon + 1 ;
    }

    if( nmuon >= 2){

      Muon *mu1 = (Muon*) branchMuon->At(0);
      Muon *mu2 = (Muon*) branchMuon->At(1);

      // Plot di-muon invariant mass 
      histDiMuonMass->Fill(((mu1->P4()) + (mu2->P4())).M());
    }
 
  }

  // Show resulting histograms
  histDiMuonMass->Write();

  f->Close();

}

void testTDime(Int_t nev = 100) {

  gSystem->Load("libEVGEN");
  gSystem->Load("libTDime");
  gSystem->Load("libdime");

  TDime* dime = new TDime();
  dime->SetEnergyCMS(7000.0);
  dime->SetYRange(-2.0, 2.0);   // Set rapidity range of mesons
  dime->SetMinPt(0.1);          // Minimum pT of mesons
  dime->Initialize();

  // (pi+pi-) histograms
  TH1* hM = new TH1D("hM", "DIME #pi^{+}#pi^{-};M_{#pi^{+}#pi^{-}} #[]{GeV/#it{c}^{2}}", 100,  0.0, 5.0);

  TClonesArray* particles = new TClonesArray("TParticle", 6);
  TParticle* part = NULL;
  TLorentzVector v[2];
  TLorentzVector vSum;

  // Event loop
  for (Int_t i = 0; i < nev; ++i) {

    dime->GenerateEvent();
    Int_t np = dime->ImportParticles(particles, "All");
    printf("\n DIME Event %d: Imported %3d particles \n", i, np);

    Int_t nPrimary = 0;

    // Loop over pion (j = 4,5) tracks
    for (Int_t j = 4; j < 6; ++j) {
      part = (TParticle*) particles->At(j); // Choose the particle
      part->Print();
      part->Momentum(v[nPrimary]);          // Copy content to v
      nPrimary++;
    }
    //particles.Clear();

    // 4-vector sum
    vSum = v[0] + v[1];

    // Fill pi+pi- histograms
    hM->Fill(vSum.M());
  }

  // Save plots as pdf
  hM->Draw();    c1->SaveAs("massTDime.pdf");

}

    void fill(const NTupleReader& tr, const double weight)
        {
            const int& nSearchBin            = tr.getVar<int>("nSearchBin");
            const double& met                = tr.getVar<double>("met");
            const double& best_had_brJet_MT2 = tr.getVar<double>("best_had_brJet_MT2");
            const int& cntCSVS               = tr.getVar<int>("cntCSVS");
            const int& nTopCandSortedCnt     = tr.getVar<int>("nTopCandSortedCnt");
            const bool& passBaseline         = tr.getVar<bool>("passBaseline");

            met_->Fill(met, weight);
            mt2_->Fill(best_had_brJet_MT2, weight);
            nt_->Fill(nTopCandSortedCnt, weight);
            nb_->Fill(cntCSVS, weight);

            if(passBaseline)
            {
                baseline_met_->Fill(met, weight);
                baseline_mt2_->Fill(best_had_brJet_MT2, weight);
                baseline_nt_->Fill(nTopCandSortedCnt, weight);
                baseline_nb_->Fill(cntCSVS, weight);

                baseline_nSearchBin_->Fill(nSearchBin, weight);
            }
        }

void TTimeHists::Fill(EHist hist)
{
   for (Long_t n = 0; n < fNum; ++n) {
      NextValues();
      if (fgDebug > 1) {
         printf("%ld: fill %s", n, hist == kHist? (fDim < 4 ? "hist" : "arr") : "sparse");
         for (Int_t d = 0; d < fDim; ++d)
            printf("[%g]", fValue[d]);
         printf("\n");
      }
      if (hist == kHist) {
         switch (fDim) {
         case 1: fHist->Fill(fValue[0]); break;
         case 2: ((TH2F*)fHist)->Fill(fValue[0], fValue[1]); break;
         case 3: ((TH3F*)fHist)->Fill(fValue[0], fValue[1], fValue[2]); break;
         default: fHn->Fill(fValue); break;
         }
      } else {
         fSparse->Fill(fValue);
      }
   }
}

void ConnectToServer(const TInetAddress *hostb, Int_t port)
{
   // Called by the Bonjour resolver with the host and port to which
   // we can connect.

   // Connect only once...
   TBonjourResolver *resolver = (TBonjourResolver*) gTQSender;
   TInetAddress host = *hostb; 
   delete resolver;

   printf("ConnectToServer: host = %s, port = %d\n", host.GetHostName(), port);

   //--- Here starts original hclient.C code ---

   // Open connection to server
   TSocket *sock = new TSocket(host.GetHostName(), port);

   // Wait till we get the start message
   char str[32];
   sock->Recv(str, 32);

   // server tells us who we are
   int idx = !strcmp(str, "go 0") ? 0 : 1;

   Float_t messlen  = 0;
   Float_t cmesslen = 0;
   if (idx == 1)
      sock->SetCompressionLevel(1);

   TH1 *hpx;
   if (idx == 0) {
      // Create the histogram
      hpx = new TH1F("hpx","This is the px distribution",100,-4,4);
      hpx->SetFillColor(48);  // set nice fillcolor
   } else {
      hpx = new TH2F("hpxpy","py vs px",40,-4,4,40,-4,4);
   }

   TMessage::EnableSchemaEvolutionForAll(gEvo);
   TMessage mess(kMESS_OBJECT);
   //TMessage mess(kMESS_OBJECT | kMESS_ACK);

   // Fill histogram randomly
   gRandom->SetSeed();
   Float_t px, py;
   const int kUPDATE = 1000;
   for (int i = 0; i < 25000; i++) {
      gRandom->Rannor(px,py);
      if (idx == 0)
         hpx->Fill(px);
      else
         hpx->Fill(px,py);
      if (i && (i%kUPDATE) == 0) {
         mess.Reset();              // re-use TMessage object
         mess.WriteObject(hpx);     // write object in message buffer
         sock->Send(mess);          // send message
         messlen  += mess.Length();
         cmesslen += mess.CompLength();
      }
   }
   sock->Send("Finished");          // tell server we are finished

   if (cmesslen > 0)
      printf("Average compression ratio: %g\n", messlen/cmesslen);

   gBenchmark->Show("hclient");

   // Close the socket
   sock->Close();
}

//.........这里部分代码省略.........
  if(fTaus.size() < 1) return true;
  cAllTauCandidates.increment();

  // Pre pT cut (for some reason in the region pT < 10 there is significantly more normal than embedded)
  for(size_t i=0; i<fTaus.size(); ++i) {
    TauCollection::Tau tau = fTaus.get(i);
    if(tau.p4().Pt() < 10) continue;
    if(!TauID::isolation(tau)) continue;
    //if(isMC() && !fIsEmbedded && tau.genMatchP4().Pt() <= 41) continue; // temporary hack
    selectedTaus.push_back(tau);
    if(!fIsEmbedded) break; // select only the first gen tau
  }
  if(selectedTaus.empty()) return true;
  cPrePtCut.increment();

  // Decay mode finding
  bool atLeastOneIsolated = false;
  for(size_t i=0; i<selectedTaus.size(); ++i) {
    TauCollection::Tau& tau = selectedTaus[i];
    if(!TauID::decayModeFinding(tau)) continue;

    if(TauID::isolation(tau)) {
      atLeastOneIsolated = true;
      hTau_AfterDecayModeFindingIsolation.fill(tau.p4(), weight);
    }
    tmp.push_back(tau);
  }
  selectedTaus.swap(tmp);
  tmp.clear();
  if(selectedTaus.empty()) return true;
  cDecayModeFinding.increment();
  if(atLeastOneIsolated) {
    if(theGenTau.isValid()) hGenTau_AfterDecayModeFindingIsolation.fill(theGenTau.p4(), weight);
    hVertexCount_AfterDecayModeFindingIsolation->Fill(fVertexCount->value(), weight);
  }

  // Eta cut
  atLeastOneIsolated = false;
  for(size_t i=0; i<selectedTaus.size(); ++i) {
    TauCollection::Tau& tau = selectedTaus[i];
    if(!TauID::eta(tau)) continue;

    if(TauID::isolation(tau)) {
      atLeastOneIsolated = true;
      hTau_AfterEtaCutIsolation.fill(tau.p4(), weight);
    }
    tmp.push_back(tau);
  }
  selectedTaus.swap(tmp);
  tmp.clear();
  if(selectedTaus.empty()) return true;
  cEtaCut.increment();
  if(atLeastOneIsolated) {
    hGenTau_AfterEtaCutIsolation.fill(theGenTau.p4(), weight);
    hVertexCount_AfterEtaCutIsolation->Fill(fVertexCount->value(), weight);
  }

  // Pt cut
  atLeastOneIsolated = false;
  for(size_t i=0; i<selectedTaus.size(); ++i) {
    TauCollection::Tau& tau = selectedTaus[i];
    if(!TauID::pt(tau)) continue;

    if(TauID::isolation(tau)) {
      atLeastOneIsolated = true;
      hTau_AfterPtCutIsolation.fill(tau.p4(), weight);

//.........这里部分代码省略.........
	if(fTree.GetLHEWeight(3) != fTree.Event_LHEWeight3) cout << "Wrong GetLHEWeight(3) Value" << endl;
	if(fTree.GetLHEWeight(4) != fTree.Event_LHEWeight4) cout << "Wrong GetLHEWeight(4) Value" << endl;
	if(fTree.GetLHEWeight(5) != fTree.Event_LHEWeight5) cout << "Wrong GetLHEWeight(5) Value" << endl;
	if(fTree.GetLHEWeight(6) != fTree.Event_LHEWeight6) cout << "Wrong GetLHEWeight(6) Value" << endl;
	if(fTree.GetLHEWeight(7) != fTree.Event_LHEWeight7) cout << "Wrong GetLHEWeight(7) Value" << endl;
	if(fTree.GetLHEWeight(8) != fTree.Event_LHEWeight8) cout << "Wrong GetLHEWeight(8) Value" << endl;
	
	// CurrentLHEWeights[0] = fTree.Event_LHEWeight /fabs(fTree.Event_LHEWeight) ;
	// CurrentLHEWeights[1] = fTree.Event_LHEWeight1/fabs(fTree.Event_LHEWeight) ;
	// CurrentLHEWeights[2] = fTree.Event_LHEWeight2/fabs(fTree.Event_LHEWeight) ;
	// CurrentLHEWeights[3] = fTree.Event_LHEWeight3/fabs(fTree.Event_LHEWeight) ;
	// CurrentLHEWeights[4] = fTree.Event_LHEWeight4/fabs(fTree.Event_LHEWeight) ;
	// CurrentLHEWeights[5] = fTree.Event_LHEWeight5/fabs(fTree.Event_LHEWeight) ;
	// CurrentLHEWeights[6] = fTree.Event_LHEWeight6/fabs(fTree.Event_LHEWeight) ;
	// CurrentLHEWeights[7] = fTree.Event_LHEWeight7/fabs(fTree.Event_LHEWeight) ;
	// CurrentLHEWeights[8] = fTree.Event_LHEWeight8/fabs(fTree.Event_LHEWeight) ;
	
	for(int i = 0 ; i<9 ;i++){
	  CurrentLHEWeights[i] = fTree.GetLHEWeight(i)/fabs(fTree.Event_LHEWeight) ;
	  WScalesTotal[i] += CurrentLHEWeights[i];
	}
      }
      #endif
      
      double cutindex = 0.5;
      bool pass=true;
      while( objcut = nextcut() ){
	ExtendedCut* thecut = (ExtendedCut*)objcut ;
	pass &= thecut->Pass(jentry , weight);
	if(! pass ){
	  break;
	}else{
	  if( isfinite (weight) ){
	    cutflowtable.Fill( cutindex , weight , EventsIsPSeudoData < fabs(weight) );
	    cutflowtablew1.Fill( cutindex , weight/fabs(weight) , EventsIsPSeudoData < fabs(weight) );
	  }
	  else
	    cout << "non-finite weight : " << weight << endl;
	}
	cutindex+=1.0;
      }

      if(! pass)
	continue;

      int cut = 0;

      //cout << alllabels[ int(cutindex) ] << endl;
//       cutflowtable.Fill( cutindex , weight , EventsIsPSeudoData < fabs(weight) );
//       cutflowtablew1.Fill( cutindex , weight/fabs(weight) , EventsIsPSeudoData < fabs(weight) );
//       cutindex ++ ;

      // if( data && !(fTree.Event_passesHLT_IsoMu20_eta2p1_v2 > 0.5) )
      // 	continue;
      
      // if( !data && !(fTree.Event_passesHLT_IsoMu20_eta2p1_v1 > 0.5) )
      // 	continue;

      nPVBeforeCutsUnCorr.Fill( fTree.Event_nPV , weight/fTree.Event_puWeight , false , true );
      nPVBeforeCuts.Fill( fTree.Event_nPV , weight , false, true);

      // if( printEventIds )
      // 	(*(EventIdFiles[cutindex])) << fTree.Event_EventNumber << endl; 
      // cutflowtable.Fill( cutindex , weight , EventsIsPSeudoData < fabs(weight) );
      // cutflowtablew1.Fill( cutindex , weight/fabs(weight) , EventsIsPSeudoData < fabs(weight) );
      // cutindex ++ ;

//.........这里部分代码省略.........
          else if(i==slimJetIdx){
            temp3Vec.SetPtEtaPhi(evt->slimJetPtVec_()[i],evt->slimJetEtaVec_()[i],evt->slimJetPhiVec_()[i]);
            NewTauJet3Vec=temp3Vec-Muon3Vec+SimTauJet3Vec;
            if(NewTauJet3Vec.Pt()>30. && fabs(NewTauJet3Vec.Eta())<2.4)HT3JetVec.push_back(NewTauJet3Vec);
            if(NewTauJet3Vec.Pt()>30. && fabs(NewTauJet3Vec.Eta())<5.)MHT3JetVec.push_back(NewTauJet3Vec);
          }

        }

        // Order the HT3JetVec and MHT3JetVec based on their pT
        HT3JetVec = utils->Order_the_Vec(HT3JetVec);
        MHT3JetVec = utils->Order_the_Vec(MHT3JetVec);


        double newHT=0,newMHT=0,newMHTPhi=-1;
        TVector3 newMHT3Vec;
        for(int i=0;i<HT3JetVec.size();i++){
          newHT+=HT3JetVec[i].Pt();
        }
        for(int i=0;i<MHT3JetVec.size();i++){
          newMHT3Vec-=MHT3JetVec[i];
        }
        newMHT=newMHT3Vec.Pt();
        newMHTPhi=newMHT3Vec.Phi();

        //New #Jet
        int newNJet = HT3JetVec.size();
        if(verbose==1)printf("newNJet: %d \n ",newNJet);



      // Acceptance determination 1: Counter for all events
      // with muons at generator level
      hAccAll->Fill( binMap_mht_nj[utils2::findBin_mht_nj(evt->nJets(),evt->mht()).c_str()] );
//      hAccAll->Fill( binMap_mht_nj[utils2::findBin_mht_nj(newNJet,newMHT).c_str()] ); // this doesn't work good

      // Check if generator-level muon is in acceptance
      if( genMuPt > LeptonAcceptance::muonPtMin() && std::abs(genMuEta) < LeptonAcceptance::muonEtaMax() ) {
        if(verbose!=0)printf("Muon is in acceptance \n ");
        // Acceptance determination 2: Counter for only those events
        // with generator-level muons inside acceptance
        // hAccPass->Fill( binMap[utils2::findBin(cntNJetsPt30Eta24,nbtag,HT,template_mht).c_str()] );
        hAccPass->Fill( binMap_mht_nj[utils2::findBin_mht_nj(evt->nJets(),evt->mht()).c_str()] );
//        hAccPass->Fill( binMap_mht_nj[utils2::findBin_mht_nj(newNJet,newMHT).c_str()] );

        // Reconstruction-efficiency determination 1: Counter for all events
        // with generator-level muons inside acceptance, regardless of whether
        // the muon has also been reconstructed or not.
        // hIsoRecoAll->Fill( binMap[utils2::findBin(cntNJetsPt30Eta24,nbtag,HT,template_mht).c_str()]);
        hIsoRecoAll->Fill( binMap[utils2::findBin_NoB(evt->nJets(),evt->ht(),evt->mht()).c_str()]);

        // Check if the muon has been reconstructed: check if a reconstructed
        // muon is present in the event that matches the generator-level muon
        // Isolation-efficiency determination 1: Counter for all events with a
        // reconstructed muon that has a generator-level muon match inside the
        // the acceptance, regardless of whether the reconstructed muon is also
        // isolated or not.

//if( evt->MuPtVec_().size()>0 )printf(" RecoMu--> Pt: %g eta: %g phi: %g deltaRMax: %g ",evt->MuPtVec_()[0],evt->MuEtaVec_()[0],evt->MuPhiVec_()[0],deltaRMax); // Ahmad3
//else cout << " Muon size is 0 \n " ;
        // in R and in pt
        int matchedMuonIdx = -1;
        if(evt->MuPtVec_().size()>0 && utils->findMatchedObject(matchedMuonIdx,genMuEta,genMuPhi,evt->MuPtVec_(),evt->MuEtaVec_(),evt->MuPhiVec_(),deltaRMax,verbose) ) {
          // Muon is reconstructed
          const double relDeltaPtMu = std::abs(genMuPt - evt->MuPtVec_().at(matchedMuonIdx) ) / evt->MuPtVec_().at(matchedMuonIdx) ;
          if(verbose!=0)printf(" relDeltaPtMu: %g \n ",relDeltaPtMu);

//.........这里部分代码省略.........
  // Get the tree from file
  TChain* tr = new TChain("AnaTree");
  tr->Add("/nfs/dust/test/cmsdas/school61/susy/ntuple/2013-v1/"+fileName(sampleId)+"_0.root");

  // Set the branches
  tr->SetBranchAddress("NrecoJet",&nRecoJets);
  tr->SetBranchAddress("recoJetPt",recoJetPt);
  tr->SetBranchAddress("recoJetPhi",recoJetPhi);
  tr->SetBranchAddress("recoJetEta",recoJetEta);
  tr->SetBranchAddress("NrecoMu",&nRecoMus);
  tr->SetBranchAddress("NrecoEle",&nRecoEle);
  if( sampleId > 0 ) tr->SetBranchAddress("EvtWgt",&evtWgt);



  // --- Process the Events in the Tree --------------------------------
  int nEvtsToProcess = tr->GetEntries();
  std::cout << "Processing " << nEvtsToProcess << " events" << std::endl;

  // Loop over the tree entries
  for(int evtIdx = 0; evtIdx < nEvtsToProcess; ++evtIdx) {
    if( evtIdx%100000 == 0 ) std::cout<<"  Event: " << evtIdx << std::endl;

    // Get the variables' values for this event
    tr->GetEntry(evtIdx);
    if( nRecoJets > kRecoJetColSize ) {
      std::cerr << "ERROR: more than " << kRecoJetColSize << " reco jets in event " << evtIdx << std::endl;
      exit(-1);
    }


    // Apply the lepton veto
    if( nRecoEle > 0 ) continue;
    if( nRecoMus > 0 ) continue;


    // Calculate RA2 selection-variables from jets
    float selNJet = 0; // Number of jets with pt > 50 GeV and |eta| < 2.5 (HT jets)
    float selHt   = 0.;
    float selMhtX = 0.;
    float selMhtY = 0.;
    // Loop over reco jets: they are ordered in pt
    for(int jetIdx = 0; jetIdx < nRecoJets; ++jetIdx) {
      // Calculate NJet and HT
      if( recoJetPt[jetIdx] > kHtJetPtMin && TMath::Abs(recoJetEta[jetIdx]) < kHtJetEtaMax ) {
	selNJet++;
	selHt += recoJetPt[jetIdx];
      }
      // Calculate MHT components
      if( recoJetPt[jetIdx] > kMhtJetPtMin && TMath::Abs(recoJetEta[jetIdx]) < kMhtJetEtaMax ) {
	selMhtX -= recoJetPt[jetIdx]*TMath::Cos(recoJetPhi[jetIdx]);
	selMhtY -= recoJetPt[jetIdx]*TMath::Sin(recoJetPhi[jetIdx]);
      }
    } // End of loop over reco jets
    float selMht = sqrt( selMhtX*selMhtX + selMhtY*selMhtY );
    

    // Select only events with at least 2 HT jets
    if( selNJet < 3 ) continue;



    //>>> PLACE OTHER RA2 CUTS HERE




    
    // Event weight in plots
    float weight = 1.;
    if( sampleId == 1 ) weight = evtWgt; // In case of the flat QCD-MC, reweight to physical spectrum

    
    // Fill histogram
    hNJets->Fill(selNJet,weight);
    hHt->Fill(selHt,weight);
    hMht->Fill(selMht,weight);
    hMEff->Fill(selHt+selMht,weight);
    for(int i = 0; i < static_cast<int>(hJetPt.size()); ++i) {
      if( i == nRecoJets ) break;
      hJetPt.at(i)->Fill(recoJetPt[i],weight);
      hJetPhi.at(i)->Fill(recoJetPhi[i],weight);
      hJetEta.at(i)->Fill(recoJetEta[i],weight);
    }
  } // End of loop over events



  // --- Save the Histograms to File -----------------------------------
  TFile outFile("General_"+fileName(sampleId)+".root","RECREATE");
  hNJets->Write();
  hHt->Write();
  hMht->Write();
  hMEff->Write();
  for(unsigned int i = 0; i < hJetPt.size(); ++i) {
    hJetPt[i]->Write();
    hJetEta[i]->Write();
    hJetPhi[i]->Write();
  }  
}

void plots() {

  gSystem->Load("libEVGEN"); // Needs to be!

  AliRunLoader* rl = AliRunLoader::Open("galice.root");
  rl->LoadKinematics();
  rl->LoadHeader();

  // 4pi histograms
  TH1* hM    = new TH1D("hM",  "DIME #rho#rho;M_{4#pi} #(){GeV/#it{c}^{2}}", 100,  1.0, 3.0);
  TH1* hPt   = new TH1D("hPt", "DIME #rho#rho;p_{T}#(){4#pi} #(){GeV/#it{c}}", 100,  0.0, 3.0);

  // pi+- histograms
  TH1* hPt1   = new TH1D("hPt1", "DIME #rho#rho;p_{T}#(){#pi^{#pm}} #(){Gev/#it{c}}",  100, 0.0, 3.0);

  AliStack* stack = NULL;
  TParticle* part = NULL;
  TLorentzVector v[4];
  TLorentzVector vSum;

  // Loop over events
  for (Int_t i = 0; i < rl->GetNumberOfEvents(); ++i) {

    rl->GetEvent(i);
    stack = rl->Stack();
    Int_t nPrimary = 0;

    // Loop over all particles
    for (Int_t j = 0; j < stack->GetNtrack(); ++j) {
      part = stack->Particle(j);           // Get particle
      part->Print();                       // Print contents

      if (abs(part->GetPdgCode()) == 211   // Is pi+ or pi-
          & part->GetStatusCode() == 1     // Is stable final state
          & stack->IsPhysicalPrimary(j)) { // Is from generator level

        part->Momentum(v[nPrimary]);       // Set content of v
        ++nPrimary;
      }
    }
    if (nPrimary != 4) {
      printf("Error: nPrimary=%d != 4 \n", nPrimary);
      continue;
    }

    // 4-vector sum
    vSum = v[0] + v[1] + v[2] + v[3];

    // Fill 4pi histograms
    hM->Fill(vSum.M());
    hPt->Fill(vSum.Perp());

    // Fill pi+- histograms
    for (Int_t k = 0; k < 4; ++k) {
      hPt1->Fill(v[k].Perp());
    }
    printf("\n");
  }

  // Save plots as pdf
  hM->Draw();    c1->SaveAs("plotM.pdf");
  hPt->Draw();   c1->SaveAs("plotPt.pdf");
  hPt1->Draw();  c1->SaveAs("plotPt1.pdf");

}