/** This method initializes histograms.
*/
INT MPulseLengths_init()
{
  // This histogram has the pulse lengths on the X-axis and the number of pulses on the Y-axis
  // One histogram is created for each detector
  
  std::map<std::string, std::string> bank_to_detector_map = gSetup->fBankToDetectorMap;
  for(std::map<std::string, std::string>::iterator mapIter = bank_to_detector_map.begin(); 
      mapIter != bank_to_detector_map.end(); mapIter++) { 

    std::string bankname = mapIter->first;
    std::string detname = gSetup->GetDetectorName(bankname);
    std::string histname = "h" + detname + "_Lengths";
    std::string histtitle = "Plot of the pulse lengths for the " + detname + " detector";
    TH1I* hDetLengths = new TH1I(histname.c_str(), histtitle.c_str(), 100,0,100);
    hDetLengths->GetXaxis()->SetTitle("Pulse Lengths [N Samples]");
    hDetLengths->GetYaxis()->SetTitle("Arbitrary Unit");
    hDetLengths->SetBit(TH1::kCanRebin);

    length_histograms_map[bankname] = hDetLengths;
  }

  std::string histname = "hAvgPulseLengthsPerChannel";
  std::string histtitle = "Plot of the average pulse lengths per event for the each channel";
  average_length_histogram = new TH2I(histname.c_str(), histtitle.c_str(), 1,0,1, 5000,0,5000);
  average_length_histogram->GetXaxis()->SetTitle("Bank Name");
  average_length_histogram->GetYaxis()->SetTitle("MIDAS Event Number");
  average_length_histogram->SetBit(TH1::kCanRebin);

  return SUCCESS;
}

bctest()
{
  // fgLegoPlotHGPtr = new AliHLTPHOSOnlineDisplayTH2D(fOnlineDisplayPtr, "Cosmics, High gain", "PHOS HLT: Cosmics", 
// 			    N_XCOLUMNS_MOD*N_MODULES , 0, N_XCOLUMNS_MOD*N_MODULES,  
// 			    N_ZROWS_MOD,   0, N_ZROWS_MOD);    

  TH1I   *plot = new TH1I("stats", "number of errors (y) vs number of transactions (x) ",20000, 0, 20001); 

  cout << "HELLO WORLD" << endl;
  FILE *fp = fopen("testresults_altro_gtl.txt", "r"); 
  int index = 0;
  int set = 0;
  int read = 0;
  int good = 0;
  int bad = 0;

  if(fp !=0)
    {
      cout << "file found" << endl; 
      
      for(int i=0; i< 20000; i++) 
	{
	  fscanf(fp,"Loop: %d\tValue set  is: 0x%xValues  read is: 0x%x,   mGood = %d,  nBad = %d\n", &index, &set, &read, &good ,&bad);
	  printf("\n index = %d,  val = %d", index, bad);
	  plot->SetBinContent(index, bad); 
	}
       
      plot->Draw();
   }
  else
    {
      cout << "Could not open file" << endl;
    }

}

void A1::PlotChannel( unsigned aModule, unsigned aChannel ){
  if (fChain == 0) return;
  if( c1 ) delete c1;
  c1 = new TCanvas( 3 );
  stringstream stmp;
  stmp << "Module " << aModule + 1 << ", Channel #" << aChannel << ";Time [0.1 ns]; Counts";


  TH1I * hChannel = new TH1I( "hChannel", stmp.str().c_str(), 100000, 10000, 100000 );
  Long64_t nentries = fChain->GetEntriesFast();
  Long64_t nbytes = 0, nb = 0;
  
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   
    nbytes += nb;
      // if (Cut(ientry) < 0) continue;
    vector<unsigned int> * HITS_PER_CHANNEL = HitsPerChannel[ aModule ];
    vector<unsigned int> * CH = _Data[ aModule ][ aChannel ];//Module2_LE_CH32;
      //std::cout << "HITS_PER_CHANNEL 32: " << HITS_PER_CHANNEL->at( 32 ) << std::endl;
    if( HITS_PER_CHANNEL->at( aChannel ) > 0 ){
      long trig = CH->at( 0 );
      hChannel->Fill( trig );
    } 
      
  }

  hChannel->Draw();

}

//______________________________________________________________________
void EventInfo(Int_t event, Int_t px, Int_t , TObject *selected)
{
   //draw the tooltip showing the backtrace for the bin at px
   if (!gTip) return;
   gTip->Hide();
   if (event == kMouseLeave)
      return;
   Double_t xpx  = gPad->AbsPixeltoX(px);
   Int_t bin = hleaks->GetXaxis()->FindBin(xpx);
   if (bin <=0 || bin > hleaks->GetXaxis()->GetNbins()) return;
   Int_t nbytes = (Int_t)hleaks->GetBinContent(bin);
   Int_t entry  = (Int_t)hentry->GetBinContent(bin);
   Int_t btid   = (Int_t)V4[entry];
   Double_t  time = 0.0001*V3[entry];
   TH1I *hbtids = (TH1I*)T->GetUserInfo()->FindObject("btids");
   if (!hbtids) return;
   if (!btidlist) btidlist = (TObjArray*)T->GetUserInfo()->FindObject("FAddrsList");
   if (!btidlist) btidlist = (TObjArray*)f->Get("FAddrsList"); //old memstat files
   if (!btidlist) return;
   Int_t nbt = (Int_t)hbtids->GetBinContent(btid-1);
   TString ttip;
   for (Int_t i=0;i<nbt;i++) {
      Int_t j = (Int_t)hbtids->GetBinContent(btid+i);
      TNamed *nm = (TNamed*)btidlist->At(j);
      if (nm==0) break;
      char *title = (char*)nm->GetTitle();
      Int_t nch = strlen(title);
      if (nch < 20) continue;
      if (nch > 100) title[100] =0;
      const char *bar = strstr(title,"| ");
      if (!bar) continue;
      if (strstr(bar,"operator new")) continue;
      if (strstr(bar,"libMemStat")) continue;
      if (strstr(bar,"G__Exception")) continue;
      ttip += TString::Format("%2d %s\n",i,bar+1);
   }
   
   if (selected) {
      TString form1 = TString::Format("  Leak number=%d, leaking %d bytes at entry=%d    time=%gseconds\n\n",bin,nbytes,entry,time);
      gTip->SetText(TString::Format("%s%s",form1.Data(),ttip.Data() ));
      gTip->SetPosition(px+15, 100);
      gTip->Reset();
   }
}

void A1_Mod::PlotCoinc( unsigned aChannelFirst, unsigned aChannelSecond ){

  if (fChain == 0) return;
  if( c1 ) delete c1;
  c1 = new TCanvas( 3 );
  stringstream stmp;
  unsigned aModule = 0;
  unsigned ChannelCoinc = 30;
  stmp << "Coincidence - Channels " << aChannelFirst << " : " << aChannelSecond << ";Time [0.1 ns]; Counts";

  
  TH1I * hChannel = new TH1I( "hChannel", stmp.str().c_str(), 10000, 0, 10000 );
  Long64_t nentries = fChain->GetEntriesFast();
  Long64_t nbytes = 0, nb = 0;
  
  for (Long64_t jentry = 0; jentry<nentries;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   
    nbytes += nb;
      // if (Cut(ientry) < 0) continue;
    vector<unsigned int> * HITS_PER_CHANNEL = HitsPerChannel[ aModule ];
    vector<unsigned int> * COINC = _Data[ aModule ][ ChannelCoinc ];
    vector<unsigned int> * TRIG = _Data[ aModule ][ 31 ];   
   //std::cout << "HITS_PER_CHANNEL 32: " << HITS_PER_CHANNEL->at( 32 ) << std::endl;
    if( (HITS_PER_CHANNEL->at( aChannelFirst ) > 0 ) && (HITS_PER_CHANNEL->at( aChannelSecond ) > 0 ) && (HITS_PER_CHANNEL->at( ChannelCoinc ) > 0 ) ){
      long ch = COINC->at( 0 );
      long trig = TRIG->at( 0 );
      std::cout << "Trig -  = " << trig - ch << std::endl;
      hChannel->Fill( trig - ch );
    } 
      
  }

  hChannel->Draw();



}

int runnestedind() {
   TTree *t = create();
   t->SetScanField(0);

   t->Scan("type[]:type[sonind[]]","type[]==3");
   TH1I* hgood = new TH1I("hgood","works",10,0,10);
   TH1I* hbad = new TH1I("hbad","used to fail",10,0,10);
   TCanvas *c = new TCanvas("c1");
   c->Divide(1,2);
   c->cd(1);
   t->Draw("type[sonind[]]>>hbad","type[]==3");
   c->cd(2);   
   t->Draw("type[sonind[]]>>hgood","type[]==3 && sonind[]>=0");

   if (hgood->GetMean()!=hbad->GetMean()) {
      cout << "Drawing \"type[sonind[]]>>hbad\",\"type[]==3\" failed!\n";
      return 1;
   }
#ifdef ClingWorkAroundErracticValuePrinter
   printf("(int)0\n");
#endif
   return 0;
}

void drawEventDifference (TDirectory* directory, TH1* refHisto, const char* fname=0)
{
  TGraphErrors* graphX = (TGraphErrors*)directory->Get("x");
  if ( graphX==0 )  return;
  TH1I* hFirst = (TH1I*)directory->Get("firstEvent");
  TH1I* hLast = (TH1I*)directory->Get("lastEvent");
  if ( hFirst==0 || hLast==0 )  return;

  std::string fullName("cDeltaE");
  if ( fname )  fullName += fname;
  else  fullName += directory->GetName();
  TCanvas* c = new TCanvas(fullName.c_str(),fullName.c_str());
  TH1* h = refHisto->Clone("DeltaE");
  h->Reset();
  h->SetTitle("DeltaE");

  TGraph* graph = new TGraph();
  graph->SetName("gDeltaE");

  double xg,yg;
  for ( unsigned int i=1; i<=hFirst->GetNbinsX(); ++i ) {
    int e1 = hFirst->GetAt(i);
    int e2 = hLast->GetAt(i);
    graphX->GetPoint(i-1,xg,yg);
    graph->SetPoint(i-1,xg,e2-e1);
  }

  double xmin,xmax,ymin,ymax;
  graph->ComputeRange(xmin,ymin,xmax,ymax);
  h->SetMinimum(0.);
  h->SetMaximum((ymax+ymin)/2.+2.*(ymax-ymin)/2.);
  h->Draw();
  graph->SetMarkerStyle(20);
//   graph->SetMarkerColor(2);
//   graph->SetLineColor(2);
  graph->Draw("P");
}

//_______________________________________________________________________
TH1I *GetNEvents(const Char_t *testfile,const Char_t *plus){
  //
  // read output
  //
  TFile *f = TFile::Open(testfile);
  if(!f || f->IsZombie()){
    printf("File not readable\n");
    return 0x0;
  }
  TString name("nEvents");
  name += plus; 
  printf("Name of nEvents %s\n",(const char*)name); 
  TH1I *l = dynamic_cast<TH1I *>(f->Get((const char*)name));
  if(!l){
    printf("nEvents not found\n");
    f->Close(); delete f;
    return 0x0;
  } 
  TH1I *returnlist = dynamic_cast<TH1I *>(l->Clone());
  if(!returnlist) return 0x0;
  returnlist->SetDirectory(0);
  f->Close(); delete f;
  return returnlist;
}

void drawTimeDifference (TDirectory* directory, TH1* refHisto, const char* fname=0)
{
  TGraphErrors* graphX = (TGraphErrors*)directory->Get("x");
  if ( graphX==0 )  return;
  TH1I* hFirst = (TH1I*)directory->Get("firstTime");
  TH1I* hLast = (TH1I*)directory->Get("lastTime");
  if ( hFirst==0 || hLast==0 )  return;

  std::string fullName("cDeltaT");
  if ( fname )  fullName += fname;
  else  fullName += directory->GetName();
  TCanvas* c = new TCanvas(fullName.c_str(),fullName.c_str());
  TH1* h = refHisto->Clone("DeltaT");
  h->Reset();
  h->SetTitle("DeltaT");

  TGraph* graph = new TGraph();
  graph->SetName("gDeltaT");

  double xg,yg;
  for ( unsigned int i=1; i<=hFirst->GetNbinsX(); ++i ) {
    std::time_t t1 = hFirst->GetAt(i);
    std::time_t t2 = hLast->GetAt(i);
    TTimeStamp ts1(hFirst->GetAt(i));
    std::cout << "Fit started at " << ts1.AsString() << std::endl;
    graphX->GetPoint(i-1,xg,yg);
    graph->SetPoint(i-1,xg,difftime(t2,t1));
  }

  double xmin,xmax,ymin,ymax;
  graph->ComputeRange(xmin,ymin,xmax,ymax);
  h->SetMinimum(0.);
  h->SetMaximum((ymax+ymin)/2.+2.*(ymax-ymin)/2.);
  h->Draw();
  graph->SetMarkerStyle(20);
//   graph->SetMarkerColor(2);
//   graph->SetLineColor(2);
  graph->Draw("P");
}

int RFMacro_BeamBunchPeriod(void)
{
	gDirectory->cd("/"); //return to file base directory

	//Goto Beam Path
	TDirectory *locDirectory = (TDirectory*)gDirectory->FindObjectAny("RF");
	if(!locDirectory)
		return 0;
	locDirectory->cd();

	//Get RF Bunch Period Histograms
	gDirectory->cd("RF_BeamBunchPeriod");
	TH1I* locHist_RFBeamBunchPeriod = (TH1I*)gDirectory->Get("RFBeamBunchPeriod");

	//Get/Make Canvas
	TCanvas *locCanvas = NULL;
	if(TVirtualPad::Pad() == NULL)
		locCanvas = new TCanvas("RF_BeamBunchPeriod", "RF_BeamBunchPeriod", 1200, 800); //for testing
	else
		locCanvas = gPad->GetCanvas();

	//Draw
	gPad->SetTicks();
	gPad->SetGrid();
	if(locHist_RFBeamBunchPeriod != NULL)
	{
		TH1I* locHist = locHist_RFBeamBunchPeriod;
		locHist->GetXaxis()->SetTitleSize(0.05);
		locHist->GetYaxis()->SetTitleSize(0.05);
		locHist->GetXaxis()->SetLabelSize(0.05);
		locHist->GetYaxis()->SetLabelSize(0.05);
		locHist->Draw();
	}

    return 1;
}

int main(int argc, char* argv[])
{
  TFile* file = new TFile("neutrons.root", "recreate");

  TChain* chain = createChain(argc, argv);
  // Assign addresses
  double target_charge, veto_charge, target_cb, veto_cb;
  unsigned long long time;
  chain->SetBranchAddress("target_total", &target_charge);
  chain->SetBranchAddress("veto_total", &veto_charge);
  chain->SetBranchAddress("target_cb", &target_cb);
  chain->SetBranchAddress("veto_cb", &veto_cb);
  chain->SetBranchAddress("time", &time);

  const int chainEntries = chain->GetEntries();
  const double adc2us = 4/1000.0;
  const int maxCharge = 30000;
  TH1F* histogram = new TH1F("timing", "timing", 2000, 0, 2000);
  TH1F* neutrons = new TH1F("nspec", "neutron spectrum", 200, 0, maxCharge);
  TH1F* hydrogen = new TH1F("nhspec", "capture on hydrogen", 200, 0, maxCharge);
  TH1F* externals = new TH1F("externspec", "externals", 200, 0, maxCharge);
  TH1F* tsa_neutrons = new TH1F("tsa_spec", "neutron spectrum from tsa plot", 200, 0, maxCharge);
  TH1I* multiplicity = new TH1I("mult", "Event multiplicity", 20, 0, 20);

  unsigned long long t1=0;
  unsigned long long t2=0;
  unsigned long long t3=0;

  // cuts
  const double max_target_cb = 0.5;
  const double max_veto_charge = 500;
  const double min_target_charge = 0;
  const double tale_min = 800;
  const double tale_max = 2000;
  const double neutrons_min = 0;
  const double neutrons_time = 50;
  const double hydro_min = 100;
  const double hydro_max = 150;

  double previous_charge = 0;

  // TSA Plots
  const int nbins = 100;
  double logmin = -0.3;
  double logmax = 4.5;
  double binwidth = (logmax-logmin)/double(nbins);
  double bin_array[nbins+1];
  bin_array[0]=std::pow(10, logmin);
  for(int i=1; i<=nbins; i++)
    bin_array[i]=bin_array[0]+std::pow(10,logmin + i*binwidth);
  TH2F* tsa = new TH2F("tsa", "time series analysis", nbins, bin_array, nbins, bin_array);

  // Keep track of timing in order to get multiplicity
  std::vector<unsigned long long> time_vec;
  time_vec.resize(1, 0);

  for(int evt=0; evt<chainEntries; ++evt)
  {
    // Print out progress
    if(!(evt%1000))
      std::cout << std::floor(double(evt)/chainEntries*100) << "%\r";
    previous_charge = target_charge;
    chain->GetEvent(evt);    

    if(veto_charge<max_veto_charge &&
       target_charge>min_target_charge &&
       target_cb < max_target_cb )
    {
      
      t3=t2;
      t2=t1;
      t1=time;
      double d12 = (t1-t2)*adc2us;
      double d23 = (t2-t3)*adc2us;
      tsa->Fill(d23, d12);
      if((time-time_vec[0])*adc2us < 50)
	time_vec.push_back(time);
      else
      {
	multiplicity->Fill(time_vec.size());
	time_vec.clear();
	time_vec.push_back(time);
      }

      histogram->Fill(d12);
      if(d12 < neutrons_time && d12 > neutrons_min)
	neutrons->Fill(target_charge);
      if(d12 < tale_max && d12 > tale_min)
	externals->Fill(target_charge);
      if(d12 < hydro_max && d12 > hydro_min)
	hydrogen->Fill(target_charge);
      if(d12 < neutrons_time && d23 < neutrons_time)
	tsa_neutrons->Fill(previous_charge);
    }
  }
  std::cout << std::endl;
  // Events in tale:
  TF1* fitter_extern = new TF1("fitter_extern", "[0]*TMath::Exp([1]*x)",
			       0, 2000);
  fitter_extern->SetParameters(1000, -4e-4);
//.........这里部分代码省略.........

int RFMacro_SelfResolution(void)
{
	gStyle->SetOptStat(1111);
	gDirectory->cd("/"); //return to file base directory

	//Goto Beam Path
	TDirectory *locDirectory = (TDirectory*)gDirectory->FindObjectAny("RF");
	if(!locDirectory)
		return 0;
	locDirectory->cd();

	//Get RF DeltaT Histograms
	gDirectory->cd("DeltaT_RF_Itself");
	TH1I* locHist_FDCRF_SelfDeltaT = (TH1I*)gDirectory->Get("FDCRF_SelfDeltaT");
	TH1I* locHist_TOFRF_SelfDeltaT = (TH1I*)gDirectory->Get("TOFRF_SelfDeltaT");
	TH1I* locHist_TAGHRF_SelfDeltaT = (TH1I*)gDirectory->Get("TAGHRF_SelfDeltaT");
	TH1I* locHist_PSCRF_SelfDeltaT = (TH1I*)gDirectory->Get("PSCRF_SelfDeltaT");

	//Time resolutions
	double locTimeResolutionSq_TOF = locHist_TOFRF_SelfDeltaT->GetStdDev() / sqrt(2.0);
	locTimeResolutionSq_TOF *= locTimeResolutionSq_TOF;
	double locTimeResolutionSq_FDC = locHist_FDCRF_SelfDeltaT->GetStdDev() / sqrt(2.0);
	locTimeResolutionSq_FDC *= locTimeResolutionSq_FDC;
	double locTimeResolutionSq_PSC = locHist_PSCRF_SelfDeltaT->GetStdDev() / sqrt(2.0);
	locTimeResolutionSq_PSC *= locTimeResolutionSq_PSC;
	double locTimeResolutionSq_TAGH = locHist_TAGHRF_SelfDeltaT->GetStdDev() / sqrt(2.0);
	locTimeResolutionSq_TAGH *= locTimeResolutionSq_TAGH;

	//Print Coarse offsets to screen:
	cout << "Time-resolution-squared for TOF/TAGH/PSC/FDC are: " << locTimeResolutionSq_TOF << ", " << locTimeResolutionSq_TAGH << ", " << locTimeResolutionSq_PSC << ", " << locTimeResolutionSq_FDC << endl;

	//Print Coarse offsets to file:
	ofstream locOutputFileStream;
	locOutputFileStream.open("rf_time_resolution_sq.txt");
	locOutputFileStream << std::setprecision(6) << locTimeResolutionSq_TOF << " " << locTimeResolutionSq_TAGH << " " << locTimeResolutionSq_PSC << " " << locTimeResolutionSq_FDC << endl;
	locOutputFileStream.close();

	//Get/Make Canvas
	TCanvas *locCanvas = NULL;
	if(TVirtualPad::Pad() == NULL)
		locCanvas = new TCanvas("RF_SelfResolution", "RF_SelfResolution", 1200, 800); //for testing
	else
		locCanvas = gPad->GetCanvas();
	locCanvas->Divide(2, 2);

	//Draw
	locCanvas->cd(1);
	gPad->SetTicks();
	gPad->SetGrid();
	if(locHist_TOFRF_SelfDeltaT != NULL)
	{
		TH1I* locHist = locHist_TOFRF_SelfDeltaT;
		locHist->GetXaxis()->SetTitleSize(0.05);
		locHist->GetYaxis()->SetTitleSize(0.05);
		locHist->GetXaxis()->SetLabelSize(0.05);
		locHist->GetYaxis()->SetLabelSize(0.05);
		locHist->Draw();
	}

	locCanvas->cd(2);
	gPad->SetTicks();
	gPad->SetGrid();
	if(locHist_TAGHRF_SelfDeltaT != NULL)
	{
		TH1I* locHist = locHist_TAGHRF_SelfDeltaT;
		locHist->GetXaxis()->SetTitleSize(0.05);
		locHist->GetYaxis()->SetTitleSize(0.05);
		locHist->GetXaxis()->SetLabelSize(0.05);
		locHist->GetYaxis()->SetLabelSize(0.05);
		locHist->Draw();
	}

	locCanvas->cd(3);
	gPad->SetTicks();
	gPad->SetGrid();
	if(locHist_PSCRF_SelfDeltaT != NULL)
	{
		TH1I* locHist = locHist_PSCRF_SelfDeltaT;
		locHist->GetXaxis()->SetTitleSize(0.05);
		locHist->GetYaxis()->SetTitleSize(0.05);
		locHist->GetXaxis()->SetLabelSize(0.05);
		locHist->GetYaxis()->SetLabelSize(0.05);
		locHist->Draw();
	}

	locCanvas->cd(4);
	gPad->SetTicks();
	gPad->SetGrid();
	if(locHist_FDCRF_SelfDeltaT != NULL)
	{
		TH1I* locHist = locHist_FDCRF_SelfDeltaT;
		locHist->GetXaxis()->SetTitleSize(0.05);
		locHist->GetYaxis()->SetTitleSize(0.05);
		locHist->GetXaxis()->SetLabelSize(0.05);
		locHist->GetYaxis()->SetLabelSize(0.05);
		locHist->Draw();
	}
}

void makehisto(TString varname, TString process, TString ofile, TString file)
{
	TFile *fo = new TFile("histos/"+ofile+".root","RECREATE");
	TFile *f = new TFile("trees_8TeV/Nominal/"+file+".root");

	TTree *tree = (TTree*)f->Get("trees/Events");

	cout << "Filling histogram " << varname << "__" << process << " with " << file << endl;

	fo->cd();

	Float_t lumi = 19739;

	// weights
	Float_t weight = 1.0;
	Float_t pu_weight = 1.0;
	Float_t btag_weight = 1.0;
	Float_t muonID_weight = 1.0;
	Float_t muonIso_weight = 1.0;
	Float_t muonTrigger_weight = 1.0;
	tree->SetBranchAddress("pu_weight", &pu_weight);
	tree->SetBranchAddress("b_weight_nominal", &btag_weight);
	tree->SetBranchAddress("muon_IDWeight", &muonID_weight);
	tree->SetBranchAddress("muon_IsoWeight", &muonIso_weight);
	tree->SetBranchAddress("muon_TriggerWeight", &muonTrigger_weight);
	// cut variables
	Float_t muonIso = 0;
	Float_t jetrms = 0;
	Float_t mtw = 0;
	tree->SetBranchAddress("mu_iso", &muonIso);
	tree->SetBranchAddress("rms_lj", &jetrms);
	tree->SetBranchAddress("mt_mu", &mtw);
	//
	Int_t wjets_class = -1;
	tree->SetBranchAddress("wjets_classification", &wjets_class);

	Float_t var = 0;
	tree->SetBranchAddress(varname, &var);
	
	TH1::SetDefaultSumw2(true);
	TH1I *hcount = (TH1I*)f->Get("trees/count_hist");
	Int_t count = hcount->GetBinContent(1); 
	cout << "sample size: " << count << endl;

	Float_t xsec = get_xsec(file);
	cout << "xsec: " << xsec << endl;
	
	// FIXME adapt range for each variable histogram
	//TH1D *histo = new TH1D(varname+"__"+process,varname+"__"+process,30,-4.5,4.5);
	//TH1D *histo = new TH1D(varname+"__"+process,varname+"__"+process,15,-4.5,4.5);
	//TH1D *histo = new TH1D(varname+"__"+process,varname+"__"+process,20,-4.5,4.5);
	// FIXME
	TH1D *histo = new TH1D(varname+"__"+process,varname+"__"+process,15,0,4.5);
	//TH1D *histo = new TH1D(varname+"__"+process,varname+"__"+process,20,0,4.5); // too much
	histo->Sumw2();
	
	Long_t nentries = tree->GetEntries();
	for(Long_t i = 0; i < nentries; i++) {
		tree->GetEntry(i);

		// wjets
		if(process == "wjets_light") {
			if(wjets_class == 0 || wjets_class == 1) continue;
		} else if (process == "wjets_heavy") {
			if(wjets_class != 0 && wjets_class != 1) continue;
		}
		
		if(jetrms >= 0.025) continue;
		if(mtw <= 50) continue;
		if(process == "qcd") {
			if (muonIso < 0.3 || muonIso > 0.5) continue;
		}
		
		if(process != "DATA"  && process != "qcd") {
			weight = pu_weight*btag_weight*muonID_weight*muonIso_weight*muonTrigger_weight;
			weight *= lumi;
			weight *= xsec;
			weight /= (Float_t)count;
		}
		
		// FIXME
		//histo->Fill(var,weight);
		histo->Fill(TMath::Abs(var),weight);
	}
	if(process == "qcd") histo->Scale(xsec);


	cout << "integral: " << histo->Integral() << endl;
	// write histos
	histo->Write();
	fo->Close();
}

//.........这里部分代码省略.........
            double adc_offset, tdc_offset;
            iss>>adc_offset>>tdc_offset; // TDC not used currently
            tagh_t_base_fadc = adc_offset;
            tagh_t_base_tdc = tdc_offset;
        }
    }
    inFile.close();

    inFile.open(inputPrefix + "fcal_base_time.txt");
    if (inFile.is_open()){
        while (getline (inFile, line)){
            istringstream iss(line);
            iss>>fcal_t_base; 
        }
    }
    inFile.close();

    inFile.open(inputPrefix + "cdc_base_time.txt");
    if (inFile.is_open()){
        while (getline (inFile, line)){
            istringstream iss(line);
            iss>>cdc_t_base; 
        }
    }
    inFile.close();


    // Do our final step in the timing alignment with tracking

    //When the RF is present we can try to simply pick out the correct beam bucket for each of the runs
    //First just a simple check to see if we have the appropriate data
    bool useRF = false;
    double RF_Period = 4.0080161;
    TH1I *testHist = Get1DHistogram("HLDetectorTiming", "TAGH_TDC_RF_Compare","Counter ID 001");
    if (testHist != NULL){ // Not great since we rely on channel 1 working, but can be craftier later.
        useRF = true;
    }
    ofstream outFile;
    TH2I *thisHist; 
    thisHist = Get2DHistogram("HLDetectorTiming", "TRACKING", "TAGM - SC Target Time");
    if (useRF) thisHist = Get2DHistogram("HLDetectorTiming", "TRACKING", "TAGM - RFBunch Time");
    if (thisHist != NULL){
        //Statistics on these histograms are really quite low we will have to rebin and do some interpolation
        outFile.open(prefix + "tagm_tdc_timing_offsets.txt", ios::out | ios::trunc);
        outFile.close(); // clear file
        outFile.open(prefix + "tagm_adc_timing_offsets.txt", ios::out | ios::trunc);
        outFile.close(); // clear file
        int nBinsX = thisHist->GetNbinsX();
        int nBinsY = thisHist->GetNbinsY();
        TH1D * selectedTAGMOffset = new TH1D("selectedTAGMOffset", "Selected TAGM Offset; Column; Offset [ns]", nBinsX, 0.5, nBinsX + 0.5);
        TH1I * TAGMOffsetDistribution = new TH1I("TAGMOffsetDistribution", "TAGM Offset; TAGM Offset [ns]; Entries", 500, -250, 250);
        for (int i = 1 ; i <= nBinsX; i++){ 
            TH1D *projY = thisHist->ProjectionY("temp", i, i);
            // Scan over the histogram
            //chose the correct number of bins based on the histogram
            float nsPerBin = (projY->GetBinCenter(projY->GetNbinsX()) - projY->GetBinCenter(1)) / projY->GetNbinsX();
            float timeWindow = 3; //ns (Full Width)
            int binWindow = int(timeWindow / nsPerBin);
            double maxEntries = 0;
            double maxMean = 0;
            for (int j = 1 ; j <= projY->GetNbinsX();j++){
                int minBin = j;
                int maxBin = (j + binWindow) <= projY->GetNbinsX() ? (j + binWindow) : projY->GetNbinsX();
                double sum = 0, nEntries = 0;
                for (int bin = minBin; bin <= maxBin; bin++){
                    sum += projY->GetBinContent(bin) * projY->GetBinCenter(bin);

void fitParamSet( string plc = "Pi", string charge = "n", int iPar = 0 ){

	gStyle->SetOptStat( 0 );
	gStyle->SetOptFit( 111 );

	XmlConfig cfg( "TpcEff.xml" );
	Logger::setGlobalLogLevel( "info" );
	vector<int> centralityBins = { 0, 1, 2, 3, 4, 5, 6 };
	vector<double> centrality = { 266.9, 218.8, 167.2, 114.3, 75.06, 36.99, 11.8 };
	vector<XmlFunction> fns;

	vector<double> x, y, ey;

	for ( int i : centralityBins ){
		INFO( "", fmt::format("{0}_{1}.TpcEffParams[{2}]", plc, charge, i) );
		XmlFunction fn( &cfg, fmt::format("{0}_{1}.TpcEffParams[{2}]", plc, charge, i) );
		fns.push_back( fn );
	}




	// Parameter A
	y.clear();
	x.clear();
	for ( int i : centralityBins ){
		y.push_back( fns[i].getTF1()->GetParameter( iPar ) );
		ey.push_back( fns[i].getTF1()->GetParError( iPar ) );
		x.push_back( centrality[ i ] );
	}
	TH1I* frame = new TH1I( "frame", (plc_label( plc, charge ) + " efficiency fit parameter " + ts( iPar ) ).c_str(), 100, 0, 300 );
	TGraphErrors * gA = new TGraphErrors( x.size(), x.data(), y.data(), 0, ey.data() );
	TF1 * fpol = new TF1( "fpol1", "pol1" );
	gA->Fit( fpol, "SE" );
	TFitResultPtr fitResult = gA->Fit( fpol, "SE" );

	double x1 = 0;
	double x2 = 300;

	TGraphErrors * band = FitConfidence::choleskyBands( fitResult, fpol, 500, 100, nullptr, x1, x2 );

	frame->Draw();
	vector<double> rangeMod = { 0.05, 0.05, 2 };
	frame->GetYaxis()->SetRangeUser( *std::min_element(y.begin(),y.end()) - rangeMod[ iPar ], *std::max_element(y.begin(),y.end()) + rangeMod[ iPar ] );
	
	gA->Draw( "same pe" );
	band->SetFillColorAlpha( kRed, 0.7 );
	band->Draw( "same e3" );
	gA->SetMarkerStyle( 8 );
	gA->SetMarkerSize( 2 );

	// gA->GetXaxis()->SetRangeUser( 0, 10 );
	// XmlFunction fn( &cfg, "Pi_n.TpcEffParams[0]" );

	
	// INFO( "", fn.getTF1()->GetParameter( 0 ) );
	

	// translate name to even format
	map<string, string> chargeMap = { { "p" , "c1" }, { "n" , "c0" } } ;
	map<string, string> plcMap = { { "Pi" , "p0" }, { "K" , "p1" }, { "P" , "p2" } };
	map<int, string> parMap = { { 0 , "a" }, { 1 , "b" }, { 2 , "c" } };


	gPad->Print( fmt::format( "efficiency_{0}_{1}_{2}.png", chargeMap[ charge ], plcMap[ plc ], parMap[ iPar ] ).c_str() );

}

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

        waveformWrap->Write();

        for(Plots p : targetPlots)
        {
            p.getTOFHisto()->Write();
            p.getEnergyHisto()->Write();
        }

        fittedTimeHisto->Write();

        deltaTVsPulseIntegral0->Write();
        deltaTVsPulseIntegral1->Write();
        deltaTVsPulseIntegral2->Write();
        deltaTVsPulseIntegral3->Write();
        deltaTVsPulseIntegral4->Write();
        deltaTVsPulseIntegral5->Write();

        deltaTVsPulseHeight->Write();

        relativeTriggerTimeHisto->Write();
        triggerAmplitudeHisto->Write();
        relativeTriggerTimeVsAmplitude->Write();

        gammaToGammaTimeH->Write();
    }

    else if(mode=="waveform")
    {
        setBranchesHistosW(treeToSort);

        triggerWalk = new TH2I("triggerWalk","trigger time vs. waveform chunk #",200,0,200,1000,0,1000);

        TH1I* monitorHisto = new TH1I("targetPosH", "targetPos", 7, 0, 7);
        monitorHisto->GetXaxis()->SetTitle("target position of each waveform");

        long triggersWithGamma = 0;
        vector<int> fullWaveform(MACRO_LENGTH/2, BASELINE);
        vector<double> triggerList;

        int prevTargetPos = 0;
        double firstTimetagInSeries = 0;

        int totalEntries = treeToSort->GetEntries();
        cout << "Total waveforms = " << totalEntries << endl;

        /*TCanvas *mycan = (TCanvas*)gROOT->FindObject("mycan");

          if(!mycan)
          {
          mycan = new TCanvas("mycan","mycan");
          }*/

        // EVENT LOOP for sorting through channel-specific waveforms
        for(int j=1; j<totalEntries; j++)
        {
            cout << "Processing triggers on waveform " << j << "\r";
            fflush(stdout);

            if(j>30)
            {
                break;
            }

            triggerList.clear();
            triggerValues.clear();

void histoTrbyTr(string decay, int saveFig= 0 , int onlyFirstCandidate= 0){
  gROOT->Reset();
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(1);
  gStyle->SetHistLineWidth(2);
  //  gStyle->SetOptStat(111110);
  //gStyle->SetOptFit(111);

  gROOT->ForceStyle();
 //  double pigreco = 3.141592;

   stringstream file;
  file.str("");
  file << "rfdir:///castor/cern.ch/cms/store/caf/user/taroni/TkAlData2011-20110721/variables"<<decay<<".root";


  cout << __LINE__ << " Opening file: " << file.str().c_str() << endl;
  
  TFile * f = new TFile();
  
  f = TFile::Open(file.str().c_str());
  if (!f){
    cout << "FILE NOT FOUND" << endl;
    return;
  }
   cout << __LINE__ << " file " << file.str().c_str() <<  " opened" << endl;
   stringstream outstring;
   outstring.str("");
   outstring << "triggerFile" << decay << ".root";

   TFile* outfile = new TFile(outstring.str().c_str(),"RECREATE");

   outfile -> cd();

//    double etaMu1;
//    double phiMu1;
//    double ptMu1;
//    double etaMu2;
//    double phiMu2;
//    double ptMu2;
   
  double nbins = 0;
  double  firstbin = 0;
  double  lastbin  = 0;

  double maxMass =0, minMass=0;

  string str3;
  size_t pos;
  pos = decay.find("MuMu");    // position of "live" in str
  str3 = decay.substr (0, pos);   // get from "live" to the end
  cout << str3 << endl;

//   if (decay == "ZMuMu" || decay == "ZMuMuMC" ) {
  if (str3 == "Z" ) {
    nbins = 140;
    firstbin = 50;
    lastbin  = 120;
    maxMass = 115;
    minMass = 65;
  }  else if (str3 == "Jpsi") {
    nbins = 80;
    firstbin = 2;
    lastbin  = 4;
    maxMass = 3.4;
    minMass = 2.7;
  }  else if (str3 == "Upsilon") {
    nbins = 80;
    firstbin = 8.5;
    lastbin  = 10.5;
    maxMass = 9.9;
    minMass = 8.9;
  }

  if (DEBUG)    cout<< __LINE__ << endl;

  TH1D * hTrigger = new TH1D("trigger","trigger",100,0.,100.); 
  
 //   TH1D *hMass= new TH1D("hMass","Invariant Mass",nbins,firstbin,lastbin);
//   TH1D *hPtMother= new TH1D("hPtMother","MotherTransverse Momentum",200,0.,200.);
//   TH1D *hEtaMother= new TH1D("hEtaMother","Mother Eta",100,-5.,5.);
//   TH1D *hPhiMother= new TH1D("hPhiMother","Mother Phi",100,-4.,4.);


  TH1D *hPt= new TH1D("hPt","Transverse Momentum",200,0.,200.);
  TH1D *hEta= new TH1D("hEta","Eta",100,-4.,4.);
  TH1D *hPhi= new TH1D("hPhi","Phi",100,-4.,4.);

//   TH1D *hDeltaPhi= new TH1D("hDeltaPhi","DeltaPhi",100,-10.,10.);
//   TH1D *hDeltaEta= new TH1D("hDeltaEta","DeltaEta",100,-10.,10.);
//   TH1D *hDeltaR  = new TH1D("R"  ,"R"  ,100,-1.,10.);

  TH1I *hMultCand = new TH1I("hMultCand","MultipleCandidate",50,-.5,49.5);

//   TH2D * hEta1vsDeltaPhi = new TH2D ("hEta1vsDeltaPhi", "EtaMu1vsDeltaPhi", 100, -2., 2. , 100, -4., 4.) ;
//   TH2D * hEta2vsDeltaPhi = new TH2D ("hEta2vsDeltaPhi", "EtaMu2vsDeltaPhi", 100, -2., 2. , 100, -4., 4.) ;

  if (DEBUG)  cout<< __LINE__ << endl;


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

void CsIProj()
{
  TFile *file = new TFile("../root/NZ_55_New.root");
  TFile *gates = new TFile("../gates/zlines.root");
  TFile *gates2 = new TFile("../gates/zlines_new.root");

  ofstream ofile("CsI_55A_New.dat");
  
  TCanvas *mycan = (TCanvas*)gROOT->FindObjectAny("mycan");
  if(!mycan)
    {
      mycan = new TCanvas("mycan","mycan");
      mycan->Divide(1,2);
    }
  
  
  ostringstream outstring;
  string name;
  int p1= 30, p2=50; //+- fit limits up to 2 peaks. May be different.
  int const num_par = 5; //number of peaks times 2(pol1)+3(gaus).
  
  for(int ic =0;ic<56;ic++)
    {
      
      outstring.str("");
      outstring << "dEE/dEE_" << ic;
      name = outstring.str();
      
      mycan->cd(1);
      TH2I *hist = (TH2I*)file->Get(name.c_str());
      hist->Draw("col");  
      hist->GetXaxis()->SetRangeUser(200.,1800.); 
      hist->GetYaxis()->SetRangeUser(5.,50.); 
   
      
      if(ic <16 || ic > 31)
	TCutG *mycut = (TCutG*)gates->Get(Form("Zline_%i_2_4",ic));
      else 
	TCutG *mycut = (TCutG*)gates2->Get(Form("Zline_%i_2_4",ic));
      mycut->Draw();

      file->cd();
      outstring.str("");
      outstring << "CsI/CsIGate/ECsI_" << ic << "_Gate";
      name = outstring.str();
      gPad->SetLogz();

      mycan->cd(2);
      TH1I * proj = (TH1I*)file->Get(name.c_str());
      proj->Draw();
      proj->Rebin(4);
      proj->GetXaxis()->SetRangeUser(700.,1800.);

      mycan->Modified();
      mycan->Update();

      TMarker * mark;
      mark=(TMarker*)mycan->WaitPrimitive("TMarker"); //Get the Background limits
      int bkg_lo = mark->GetX();
      delete mark;  
      mark=(TMarker*)mycan->WaitPrimitive("TMarker");
      int bkg_hi = mark->GetX();
      delete mark;
      mark=(TMarker*)mycan->WaitPrimitive("TMarker"); // Get the 1st peak initial guess
      int peak1 = mark->GetX();
      delete mark;
      
      
      double par[num_par] = {0.};
      double out[num_par] = {0.}; 
      int peak1_lo = peak1 - p1, peak1_hi = peak1 + p1; // Peak center and limits
      
      
      TF1 *l1 = new TF1("l1", "pol1", bkg_lo, bkg_hi);
      TF1 *g1 = new TF1("g1", "gaus", peak1_lo,peak1_hi);
      
      TF1 *total = new TF1("total", "pol1(0)+gaus(2)", bkg_lo,bkg_hi);
      
      proj->Fit(l1,"R");
      proj->Fit(g1,"R+");
      
      l1->GetParameters(&par[0]);
      g1->GetParameters(&par[2]);
      
      total->SetParameters(par);
      proj->Fit(total,"R");
      total->GetParameters(out);
      
      
      ofile << ic << " " << out[3] << endl;
      
      outstring.str("");
      outstring << "55A_" << ic;
      name = outstring.str();
      total->SetName(name.c_str());
      total->Draw("same");
      mycan->Modified();
      mycan->Update();
      
      bool IsGood = 0;
//.........这里部分代码省略.........