C++ TH1F::GetMaximumBin方法代码示例

bool isgood(TH1F h, float sigma){
  if(sigma == -1) return true;
  bool Flag = true;
  int imaxbin = h.GetMaximumBin();
  float bincenter = h.GetBinCenter(imaxbin);
  for(int ibin=0;ibin<h.GetNbinsX();ibin++){
      if((h.GetBinCenter(ibin)<(bincenter-sigma)) || (h.GetBinCenter(ibin)>(bincenter+sigma))){
          if(h.GetBinContent(ibin)!=0){
              Flag = false;
              break;
          }
      }
  }

  if(Flag){
  h.SetBinContent(imaxbin,0);
  imaxbin = h.GetMaximumBin();
  bincenter = h.GetBinCenter(imaxbin);
  for(int ibin=0;ibin<h.GetNbinsX();ibin++){
      if((h.GetBinCenter(ibin)<(bincenter-sigma)) || (h.GetBinCenter(ibin)>(bincenter+sigma))){
          if(h.GetBinContent(ibin)!=0){
              Flag = false;
              break;
          }
      }
  }
  }
  return Flag;
}

void VetoTimes()
{
  TH1F* Time;
  Char_t Buff[256];

  for(Int_t ch=0; ch<384; ch++)
  {
    sprintf(Buff, "Veto_Time%d", ch);
    Time = (TH1F*)gROOT->FindObject(Buff);
    printf("%6.1f\n", (Time->GetMaximumBin()*Time->GetBinWidth(Time->GetMaximumBin()) - 500.0)/0.05000);
  }
}

void PIDPedes()
{
  TH1F* Energy;
  Char_t Buff[256];

  for(Int_t ch=0; ch<24; ch++)
  {
    sprintf(Buff, "PID_Energy%d", ch);
    Energy = (TH1F*)gROOT->FindObject(Buff);
    printf("%4.0f\n", Energy->GetMaximumBin()*Energy->GetBinWidth(Energy->GetMaximumBin()));
  }
}

void PIDTimes()
{
  TH1F* Time;
  Char_t Buff[256];

  for(Int_t ch=0; ch<24; ch++)
  {
    sprintf(Buff, "PID_Time%d", ch);
    Time = (TH1F*)gROOT->FindObject(Buff);
    printf("%5.1f\n", Time->GetMaximumBin()*Time->GetBinWidth(Time->GetMaximumBin()) - 200.0);
  }
}

void LadderTimes()
{
  TH1F* Time;
  Double_t par[4];
  Double_t Max;
  Char_t Buff[256];

  for(Int_t ch=0; ch<352; ch++)
  {
    sprintf(Buff, "Ladder_Time%d", ch);
    Time = (TH1F*)gROOT->FindObject(Buff);
    Max = Time->GetMaximumBin()*Time->GetBinWidth(Time->GetMaximumBin()) - 1000.0;
   
    TF1* gauss = new TF1("gauss", "gaus", Max-8, Max+8);
    Time->Fit(gauss, "RQ+");
    gauss->GetParameters(&par[0]);
    
    printf("%5.2f\n", par[1]/0.117710);
  }
}

void AlignTopAndBottomTOFs(string ArrayPart="NV",Float_t Start=-5,Float_t End=25){
  
  for (int bar=1;bar<13;bar++){
    stringstream nameTop;
    stringstream nameBottom;

    nameTop<<ArrayPart<<setfill('0')<<setw(2)<<bar<<"_TopCutTOFPr";
    TH1F * CurrentTop = (TH1F*)gDirectory->Get(nameTop.str().c_str());


    nameBottom<<ArrayPart<<setfill('0')<<setw(2)<<bar<<"_BottomCutTOFPr";
    TH1F * CurrentBottom = (TH1F*)gDirectory->Get(nameBottom.str().c_str());
  
    // TCanvas *c=new TCanvas("c");
    // c->cd(1);
    // CurrentTop->Draw();
    // CurrentBottom->Draw("same");

    int binsTop = CurrentTop->GetNbinsX();

    int binsBottom= CurrentBottom->GetNbinsX();
  
    if ( binsTop !=binsBottom){
      cout<<"Nubmer of bins for the top channel not the same as for the bottom channel"<<endl;
      cout<<"This was for "<<CurrentTop->GetName()<<" and "<<CurrentBottom->GetName()<<endl;
      return;
    }
  
    int bins=binsTop;
  
    double topLow = CurrentTop->GetBinLowEdge(1);//ROOT histograms start at bin 1. bin 0 is underflow
    double topHigh= CurrentBottom->GetBinLowEdge(binsTop) +CurrentBottom->GetBinWidth(binsTop);//ROOT histograms end at bin N. N+1 is overflow bin.

    double NanoSecsPerBin = ((topHigh-topLow)*4)/binsTop;
  
    //    cout<<topLow<<" "<<topHigh<<endl;
    //    cout<<"NanoSecsPerBin "<<NanoSecsPerBin<<endl;
  
 
    int ForthOfBins = TMath::Floor(0.25* bins);
    int EighthOfBins = TMath::Floor( (1.0/8)* bins);

    int NumberOfShifts = 2*EighthOfBins; // 1/8 of the number of bins *2 for left/right shifts
  
    int TopZeroBin = CurrentTop->GetMaximumBin();
  
    int StartBin = CurrentTop->FindBin(Start);
    int EndBin = CurrentTop->FindBin(End);

    //Find shift That moves Bottom -> Top
    vector <double> TheChi2s;
    vector <int> TheShifts;
    for (int i=0;i<NumberOfShifts;i++){
      double binShift = (i - NumberOfShifts/2);

      double chi2=0;
      // cout<<"Looking in bin range "<<TopZeroBin-EighthOfBins<<" "<<TopZeroBin+EighthOfBins<<endl;
      // cout<<"That is from "<<CurrentTop->GetBinCenter(TopZeroBin-EighthOfBins)<<" "<<CurrentTop->GetBinCenter(TopZeroBin+EighthOfBins)<<endl;

      for (int bin=StartBin ;bin<EndBin;bin++){
	double b = CurrentBottom->GetBinContent(bin + binShift);
	double t = CurrentTop->GetBinContent(bin);
	if (b !=0 && t!=0){
	  double temp =((b-t)*(b-t))/( TMath::Sqrt(t) );
	  chi2+=temp;
	}
      }
    
      TheChi2s.push_back(chi2);
      TheShifts.push_back(binShift);

    }
  
    ///Find miminum by linear search 
    double min = 9999999999.0;
    int MinSpot=-1;
    for (int i=0;i<NumberOfShifts;i++){
      if ( TheChi2s[i] < min){
	min=TheChi2s[i];
	MinSpot=i;
      }
    }

    cout<<ArrayPart<<setfill('0')<<setw(2)<<bar<<"B  "<<"1  0  "<<fixed<<setw(6)<<setprecision(4)<<TheShifts[MinSpot]/4.0*NanoSecsPerBin<<endl;
    cout<<ArrayPart<<setfill('0')<<setw(2)<<bar<<"T  "<<"1  0  0"<<endl;

  }


  cout<<"There were "<<NanoSecsPerBin<<" nano secs per bin"<<endl;
  // TGraph * graph = new TGraph();
  
  // cout<<"Size of Chi2s "<<TheChi2s.size()<<endl;
  // cout<<"Num "<<NumberOfShifts<<endl;

  // for (int i=0;i<NumberOfShifts;i++){
  //   graph->SetPoint(i,TheShifts[i],TheChi2s[i]);
  // }

  // TCanvas *c2 = new TCanvas("c2");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	hToTMap->DrawCopy("colz");

  
  std::cout<<"Making per pixel tot map"<<std::endl;
  // Per pixel map
  TH2F* hToTMapPerPixel = new TH2F("hToTMapPerPixel","hToTMapPerPixel",nPixelsX,0,nPixelsX,nPixelsY,0,nPixelsY);
  for(int x=0;x<nPixelsX;x++){
    for(int y=0;y<nPixelsY;y++){
      int id = x + y*nPixelsX;
      // Get the ToT distribution for this location
      string histoname = "/tpanal/ClicpixAnalysis/hClusterTOTAssociated1pixPixel"+makestring(id);
      TH1F* totMapHisto = 0;
      totMapHisto = (TH1F*)gDirectory->Get(histoname.c_str());
      if(totMapHisto == 0) continue;
      if(totMapHisto->GetEntries() < 10) continue;
      
      // Fit the tot distribution and enter the value in the ToT map
      totMapHisto->Fit("gaus","0q");
      TF1* fit = totMapHisto->GetFunction("gaus");
      hToTMapPerPixel->Fill(x,y,fit->GetParameter(1));
      
      delete fit; delete totMapHisto;
    }
  }
  // Plot the map of ToT
  TCanvas* canv2 = new TCanvas();
  hToTMapPerPixel->GetXaxis()->SetTitle("Column");
  hToTMapPerPixel->GetYaxis()->SetTitle("Row");
  hToTMapPerPixel->GetZaxis()->SetTitle("Mean charge (ToT)");
  hToTMapPerPixel->SetMaximum(10);
  hToTMapPerPixel->SetMinimum(0);
  hToTMapPerPixel->DrawCopy("colz");

  std::cout<<"Making pixel response tot maps"<<std::endl;
  // Pixel response maps
  int responseBins=100;
  double responseWidth=0.05;
  TH2F* hToTMapPixelResponseX = new TH2F("hToTMapPixelResponseX","hToTMapPixelResponseX",responseBins+1,-responseWidth,responseWidth,15,0,15);
  for(double rID=0;rID<responseBins;rID++){
    // Get the ToT distribution for this location
    string histoname = "/tpanal/ClicpixAnalysis/hClusterTOTAssociated1pixResponseX"+makestringD(rID);
    TH1F* totMapHisto = 0;
    totMapHisto = (TH1F*)gDirectory->Get(histoname.c_str());
    if(totMapHisto == 0) continue;
    if(totMapHisto->GetEntries() < 1) continue;
 		// loop over all tot bins
    for(int i=0;i<45;i++){
      hToTMapPixelResponseX->Fill(2.*responseWidth*(rID/100.)-responseWidth,i,totMapHisto->GetBinContent(totMapHisto->FindBin(i)));
    }
    delete totMapHisto;
  }
  // Plot the map of ToT
  TCanvas* canv3 = new TCanvas();
  hToTMapPixelResponseX->GetXaxis()->SetTitle("Track distance x (mm)");
  hToTMapPixelResponseX->GetYaxis()->SetTitle("ToT");
  hToTMapPixelResponseX->DrawCopy("colz");

  TH1F* hToTProfilePixelResponseY = new TH1F("hToTProfilePixelResponseY","hToTProfilePixelResponseY",responseBins+1,-responseWidth,responseWidth);
  TH2F* hToTMapPixelResponseY = new TH2F("hToTMapPixelResponseY","hToTMapPixelResponseY",responseBins+1,-responseWidth,responseWidth,16,0,16);
  for(double rID=0;rID<responseBins;rID++){
    // Get the ToT distribution for this location
    string histoname = "/tpanal/ClicpixAnalysis/hClusterTOTAssociated1pixResponseY"+makestringD(rID);
    TH1F* totMapHisto = 0;
    totMapHisto = (TH1F*)gDirectory->Get(histoname.c_str());
    if(totMapHisto == 0) continue;
    if(totMapHisto->GetEntries() < 1) continue;
    // loop over all tot bins
    for(int i=0;i<45;i++){
      hToTMapPixelResponseY->Fill(2.*responseWidth*(rID/100.)-responseWidth,i,totMapHisto->GetBinContent(totMapHisto->FindBin(i)));
    }
    hToTProfilePixelResponseY->Fill(2.*responseWidth*(rID/100.)-responseWidth+0.001,totMapHisto->GetXaxis()->GetBinCenter(totMapHisto->GetMaximumBin()));
    delete totMapHisto;
  }
  // Plot the map of ToT
  TCanvas* canv4 = new TCanvas();
  hToTMapPixelResponseY->GetXaxis()->SetTitle("Track distance y (mm)");
  hToTMapPixelResponseY->GetYaxis()->SetTitle("ToT");
  hToTMapPixelResponseY->GetZaxis()->SetTitle("Number of events");
  hToTMapPixelResponseY->DrawCopy("colz");
//  drawline(-0.03,1,7);
//  drawline(0.03,1,7);
	
  TCanvas* canv6 = new TCanvas();
  hToTProfilePixelResponseY->GetXaxis()->SetTitle("Track distance y (mm)");
  hToTProfilePixelResponseY->GetYaxis()->SetTitle("Most probable ToT");
  hToTProfilePixelResponseY->DrawCopy("");



  
  TCanvas* canv5 = new TCanvas();
  TH1F* sidebandLow = (TH1F*)gDirectory->Get("/tpanal/ClicpixAnalysis/hClusterTOTAssociated1pixResponseY20");
  TH1F* sidebandHigh = (TH1F*)gDirectory->Get("/tpanal/ClicpixAnalysis/hClusterTOTAssociated1pixResponseY80");
  sidebandHigh->GetXaxis()->SetTitle("ToT (single pixel)");
  sidebandHigh->GetXaxis()->SetRangeUser(0,20);
  sidebandHigh->DrawCopy("");
  sidebandLow->SetLineColor(2);
  sidebandLow->DrawCopy("same");

}

//.........这里部分代码省略.........
  hQCD->Scale(kfactor);

  TH1F *hBkg = (TH1F*)hQCD->Clone("hBkg");
  hBkg->Add(h[1]);
  hBkg->Add(h[2]);
  hBkg->Add(h[3]);
  hBkg->Add(h[4]);
  hBkg->Add(hST);
  //hBkg->SetFillColor(kGray);

  cout<<"======== "<<VAR<<"====================="<<endl;
  cout<<"Data events:  "<<h[0]->Integral()<<endl;
  cout<<"QCD events:   "<<hQCD->Integral()<<endl;
  cout<<"WJets events: "<<h[2]->Integral()<<endl;
  cout<<"ZJets events: "<<h[3]->Integral()<<endl;
  cout<<"ST events:    "<<hST->Integral()<<endl;
  cout<<"TTbar events: "<<h[1]->Integral()<<endl;
  cout<<"kfactor:      "<<kfactor<<endl;

  THStack *hs = new THStack("hs","hs");
  if (LOG) {   
    hs->Add(h[2]);
    hs->Add(h[3]);
    hs->Add(hST); 
    hs->Add(hQCD);
    hs->Add(h[1]);
  }
  else {
    hs->Add(h[3]);
    hs->Add(hST);
    hs->Add(h[2]);
    hs->Add(hQCD);
    hs->Add(h[1]);
  }
  TH1F *hRatio = (TH1F*)h[0]->Clone("Ratio");
  hRatio->SetLineWidth(2);
  hRatio->Divide(hBkg);

  TLegend *leg = new TLegend(0.86,0.7,0.99,0.9);
  leg->SetFillColor(0);
  leg->SetTextFont(42);
  leg->SetTextSize(0.03);
  leg->AddEntry(hQCD,"QCD","F");
  leg->AddEntry(h[1],"TTbar","F"); 
  leg->AddEntry(hST,"ST","F"); 
  leg->AddEntry(h[2],"WJets","F");
  leg->AddEntry(h[3],"ZJets","F"); 
  
  can->SetBottomMargin(0.25);
  TH1F *hAux = (TH1F*)h[0]->Clone("aux");
  hAux->Reset();
  hAux->GetXaxis()->SetNdivisions(XNDIV);
  if (isINT) {
    hAux->GetXaxis()->CenterLabels();
  } 
  hAux->GetYaxis()->SetRangeUser(0.5,1.1*TMath::Max(hBkg->GetBinContent(hBkg->GetMaximumBin()),h[0]->GetBinContent(h[0]->GetMaximumBin())));  
  if (LOG) {
    gPad->SetLogy(); 
    hAux->GetYaxis()->SetRangeUser(0.5,2*TMath::Max(hBkg->GetBinContent(hBkg->GetMaximumBin()),h[0]->GetBinContent(h[0]->GetMaximumBin())));
  }
  hAux->GetXaxis()->SetRangeUser(XMIN,XMAX);
  hAux->GetYaxis()->SetTitle(TString::Format("Number of events / %1.2f fb^{-1}",LUMI/1000));
  hAux->GetXaxis()->SetTitle("");
  hAux->GetXaxis()->SetLabelSize(0.0);
  hAux->Draw();
  hs->Draw("hist same");
  //hBkg->Draw("sames hist"); 
  h[0]->Draw("sames E");
  
  leg->Draw();
  gPad->RedrawAxis();

  TPad *pad = new TPad("pad","pad",0.,0.,1.,1.);
  pad->SetTopMargin(0.77);
  pad->SetRightMargin(0.15);
  pad->SetFillColor(0);
  pad->SetFillStyle(0);
  pad->Draw();
  pad->cd(0);
  pad->SetGridy();
  hRatio->GetXaxis()->SetTitleOffset(0.95);
  hRatio->GetYaxis()->SetTitleOffset(1.5);
  hRatio->GetYaxis()->SetTickLength(0.06);
  hRatio->GetYaxis()->SetTitleSize(0.03);
  hRatio->GetYaxis()->SetLabelSize(0.03);
  hRatio->GetYaxis()->SetTitle("Data/MC");
  hRatio->GetXaxis()->SetTitle(XTITLE);
  hRatio->GetXaxis()->SetRangeUser(XMIN,XMAX);
  hRatio->GetYaxis()->SetRangeUser(0.5,1.5);
  hRatio->GetYaxis()->SetNdivisions(505);
  hRatio->GetXaxis()->SetNdivisions(XNDIV);
  if (isINT) {
    hRatio->GetXaxis()->CenterLabels();
  }
  hRatio->Draw();
  if (PRINT) {
    can->Print("plots/"+TString(can->GetName())+".pdf"); 
    can->Print("plots/"+TString(can->GetName())+".png");
  }
}

void checkflat(int icent = 0, int ihar=0){
    int color[nsub]={1,2,3,4,5,6,7,8};
    float pi = acos(-1.0);
    TString str;
    TFile *fin;
        TH1::SetDefaultSumw2();
        gStyle->SetOptStat(kFALSE);
        int iharE=0;
     if(nhar==1||nhar==2) iharE=1;
        int n = ihar+1.0+iharE;
    int nrun = GetTotalRun();
    float FitGood[nsub][1000];
    float FitGooderr[nsub][1000];
    float runlist[1000];
        TH1F* hpsi = new TH1F("psi","psi",100,-pi,pi);
        TH1F* hunpsi = new TH1F("unpsi","unpsi",100,-pi,pi);
        TH1F* h = new TH1F("h","",100,-4,4);
        for(int irun=0;irun<nrun;irun++){
        std::ifstream corrs("run15pAl200MBPro104.Lst");
        int index=0; int run=0;
        for(int jrun=0;jrun<irun+1;jrun++){
        corrs>>index>>run;
        }
        runlist[irun] = run;
       // system(Form("hadd -f /phenix/plhf/xuq/taxi/Run15pAl200MBPro104/10241/data/%d.root /phenix/plhf/xuq/taxi/Run15pAl200MBPro104/10241/data/%d_*.root",run,run));
        fin = TFile::Open(Form("/phenix/plhf/xuq/taxi/Run15pAl200MBPro104/10241/data/%d.root",run));
        for(int isub=0;isub<nsub;isub++){
        str = choosesub(isub);
        if(str=="ABORT") continue;
        hpsi->Reset();
        hunpsi->Reset();
        for(int ibbcz=0;ibbcz<nbbcz;ibbcz++){
          TH1F* hunpsitemp = (TH1F*)fin->Get(Form("psi_0_0_%d_%d_%d_%d",icent,ibbcz,ihar,isub));
          TH1F* hpsitemp = (TH1F*)fin->Get(Form("psiFla_0_0_%d_%d_%d_%d",icent,ibbcz,ihar,isub));
          hpsi->Add(hpsitemp);
          hunpsi->Add(hunpsitemp);
        }
        TF1 *fun = new TF1("fun","pol0",-pi,pi);
        TCanvas *c1 = new TCanvas("c1","c1",600,450);
        h->GetXaxis()->SetTitle(Form("#Psi_{%d}",n));
        h->GetYaxis()->SetTitle(Form("# of events"));
        h->GetYaxis()->SetRangeUser(hunpsi->GetBinContent(hunpsi->GetMinimumBin())*0.9,hunpsi->GetBinContent(hunpsi->GetMaximumBin())*1.1);
        h->Draw("C");
       // h->GetYaxis()->SetRangeUser(0,hunpsi->GetBinContent(hunpsi->GetMaximumBin())*1.2);
      if(hpsi->GetEntries()>1000){
	hpsi->Fit("fun","QR0");
	float par=fun->GetParameter(0);
	float parerr=fun->GetParError(0);
        FitGood[isub][irun]=fun->GetChisquare()/fun->GetNDF();
        FitGooderr[isub][irun]=0;
        fun->SetLineColor(4);
        fun->Draw("same");
      }
        else{
        FitGood[isub][irun]=0;
        FitGooderr[isub][irun]=0;
        }
        if(irun!=104) continue;
	hpsi->SetMarkerStyle(20);
	hpsi->SetMarkerSize(0.8);
	hpsi->SetMarkerColor(2);
	hpsi->SetLineColor(2);
	hunpsi->SetMarkerStyle(20);
	hunpsi->SetMarkerSize(0.8);
	hunpsi->SetMarkerColor(1);
	hunpsi->SetLineColor(1);
        TLegend *leg = new TLegend(0.4,0.7,0.8,0.85);
        leg->SetTextSize(0.04);
        leg->SetFillColor(0);
        leg->SetBorderSize(0);
        leg->AddEntry(hunpsi,Form("Run%d %s",run,str.Data()));
        leg->AddEntry(hunpsi,"before flattening","L");
        leg->AddEntry(hpsi,"After flattening","L");
        hpsi->Draw("Csame");
        hunpsi->Draw("Csame");
        leg->Draw("same");
        c1->Print(Form("run-by-run/checkflat_%d_%d_%d_%d.png",icent,ihar,isub,irun));
        delete c1;
        }
        fin->Close();
    }  
        TGraphErrors *grRun[nsub];
    TCanvas *c1 = new TCanvas("c1","c1",600,450);
    TCanvas *c2 = new TCanvas("c2","c2",600,450);
    TH1F* h1 = new TH1F("h1","",4000,436000,439000);
    h1->GetXaxis()->SetTitle(Form("Run Number"));
    h1->GetYaxis()->SetTitle("Fit #chi^{2}/Ndf");
    h1->GetXaxis()->SetRangeUser(runlist[0]-15,runlist[nrun-1]+15);
    h1->GetYaxis()->SetRangeUser(0,4);
    TLegend *leg = new TLegend(0.65,0.75,0.88,0.85);
    leg->SetTextSize(0.035);
    leg->SetBorderSize(0);
    leg->SetFillColor(0);
    c1->cd();
    h1->Draw("0");
        for(int isub=0;isub<nsub;isub++){
            if(isub!=4 && isub!=5) continue;
            grRun[isub] = new TGraphErrors(nrun,runlist,FitGood[isub],0,FitGooderr[isub]);
            grRun[isub] -> SetMarkerColor(isub-3);
            grRun[isub] -> SetMarkerSize(1.);
//.........这里部分代码省略.........

void calCrystalDepo(TTree *t1041, int ThisModule = -99, bool IsBatch =0){
  if(IsBatch) gROOT->SetBatch();

  TString ThisModuleString = TString::Itoa(ThisModule,10);
  int UpOrDown = ThisModule/abs(ThisModule);
  float Scale = 1.;
  if(UpOrDown == -1) Scale = 0.5;
  cout << "Creating energy distributions of events with maximum deposition in module: " << ThisModuleString << endl;
//Define Histograms

//This histogram is a tool
  TH1F* MatrixDepo = new TH1F("MatrixDepo", "Deposition on Matrix", 68, -17, 17);

  TH1F* fracDep = new TH1F("fracDep", "Highest Amplitude Crystal/Total Deposition", 100, 0, 1.2);
  TH1F* MaxBin = new TH1F("MaxBin", "Deposition of Highest Amplitude Crystal", 100,400*Scale,7000*Scale);
  TH1F* Integral = new TH1F("Integral", "Integral", 100,4000*Scale,13000*Scale);
  TH1F* SigNoise = new TH1F("SigNoise", "Signa/Pedestal", 100, 0, 2.);

//

  gStyle->SetOptStat(0);
  

  Mapper *mapper=Mapper::Instance();

  TBEvent *event = new TBEvent();
  t1041->SetBranchAddress("tbevent", &event);
  bool haverechits=false;
  vector<TBRecHit> *rechits=0;

  if(t1041->GetListOfBranches()->FindObject("tbrechits")) {
    cout <<"found rechits"<<endl;
    t1041->SetBranchAddress("tbrechits",&rechits);
    haverechits=true;
  }

  Int_t start=0; Int_t end=t1041->GetEntries();  
  
  int nEntries=0;
//  end = 10;
  for (Int_t i=start; i<end; i++) {
    t1041->GetEntry(i);
    if (i==0) mapper->SetEpoch(event->GetTimeStamp());
    if(haverechits && rechits->size() < 1) continue;

    MatrixDepo->Reset();
    float maxDepoModuleID = -99;
    float maxDepo = 0;
    float secondMaxDepo = 0;
    float secondMaxDepoModuleID = -99;
    float totdep = 0;
    float totalnoise = 0;
    float sumrms = 0;
    for (Int_t j = 0; j < event->NPadeChan(); j++){
      if (haverechits && j>=(int)rechits->size()) break;
      double ped,sig, max, maxTime;
      int channelID;
      if (haverechits){
	TBRecHit &hit=rechits->at(j);
	ped=hit.Pedestal();
	sig=hit.NoiseRMS();
	max=hit.AMax();
	maxTime=hit.TRise();
	channelID=hit.GetChannelID();
      }

      double x,y;
      int moduleID,fiberID;
      mapper->ChannelID2ModuleFiber(channelID,moduleID,fiberID);
      mapper->ModuleXY(moduleID,x,y);

      MatrixDepo->Fill(moduleID, max);      
      if(UpOrDown > 0 && moduleID > 0)	totalnoise+=ped;
      if(UpOrDown < 0 && moduleID < 0)	totalnoise+=ped;
    }

    int minx = 0;
    int maxx = 0;
    if(UpOrDown > 0){
      minx = MatrixDepo->GetXaxis()->FindBin(0.);
      maxx = MatrixDepo->GetXaxis()->FindBin(17.);
      MatrixDepo->GetXaxis()->SetRange(minx, maxx);
    }
    if(UpOrDown < 0){
      minx = MatrixDepo->GetXaxis()->FindBin(-17.);
      maxx = MatrixDepo->GetXaxis()->FindBin(0.);	
      MatrixDepo->GetXaxis()->SetRange(minx, maxx);
    }

//    cout << "TOTAL NOISE: " << totalnoise << endl;

    float maxbin = MatrixDepo->GetBinContent(MatrixDepo->GetMaximumBin());
    float integral = MatrixDepo->Integral(minx,maxx);

    maxDepoModuleID = MatrixDepo->GetBinLowEdge(MatrixDepo->GetMaximumBin());

    if( abs(ThisModule) != 99 && maxDepoModuleID != ThisModule ) continue;
    if( abs(ThisModule) == 100 ){
	if( abs(maxDepoModuleID) != 6 &&
	    abs(maxDepoModuleID) != 7 &&
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
  if (argc < 4)
    {
      cout << "Usage: ./simData [numfiles] [Elow] [Ehigh] [geometry]\n";
      return -1;
    }

  Char_t temp2[500];
  ofstream outfile0, outfile1;
  
  Int_t num_files = atoi(argv[1]);
  Int_t Emin = atoi(argv[2]);
  Int_t Emax = atoi(argv[3]);
  Char_t *geometry = argv[4];  

  //Int_t EnergyRange = 800;
  //Int_t nbins = EnergyRange/binWidthEnergy;
  Int_t nbins = 1; //this is just remnant for now.
  //Double_t Energy_min = 0., Energy_max = (double)binWidthEnergy;
  Char_t temp[200], temp1[100], file0[100], file1[100];
  Char_t *calDir = getenv("UCNA_CAL_DIR"); 
  Char_t *G4Dir = getenv("G4WORKDIR");
  sprintf(file0, "%s/Sims/Beta_%i-%ikev_%s_east.dat", calDir, Emin, Emax, geometry);
  sprintf(file1, "%s/Sims/Beta_%i-%ikev_%s_west.dat", calDir, Emin, Emax, geometry);
  cout << file0 << endl;
  outfile0.open(file0);
  outfile1.open(file1);


  TChain *chain = new TChain("anaTree");

  
      
  for (Int_t file_number = 0; file_number<num_files; file_number++)
    {
	  
      sprintf(temp2, "%s/output/%sgeo_n1_1mil/analyzed_%i.root", G4Dir, geometry, file_number);
      chain->Add(temp2);
    }


  TH1F *east = new TH1F("east", "East Side", 200, 0., 50.);
  TH1F *west = new TH1F("west", "West Side", 200, 0., 50.); 

  for (int i = 0; i<nbins; i++)
    {
      east->Reset();
      west->Reset();
	  

      sprintf(temp, "MWPCEnergy.MWPCEnergyE>1e-5 && MWPCEnergy.MWPCEnergyE<50. && MWPCEnergy.MWPCEnergyW<0.01 && Edep.EdepW<0.1 && Edep.EdepE>%f && Edep.EdepE<%f",(double)Emin,(double)Emax);
      sprintf(temp1, "MWPCEnergy.MWPCEnergyE>>east");
      chain->Draw(temp1, temp);
	  

      sprintf(temp, "MWPCEnergy.MWPCEnergyW>1e-5 && MWPCEnergy.MWPCEnergyW<50. && MWPCEnergy.MWPCEnergyE<0.01 && Edep.EdepE<0.1 && Edep.EdepW>%f && Edep.EdepW<%f",(double)Emin,(double)Emax);
      sprintf(temp1, "MWPCEnergy.MWPCEnergyW>>west");
      chain->Draw(temp1, temp);

	 
      TF1 *f1 = new TF1("f1", "landau", 0.05, 20.);
      //f1->SetParLimits(0, 0., 1.e10);
      //f1->SetParLimits(1, 0., 1.e7);
      
      Double_t max_bin = east->GetBinCenter(east->GetMaximumBin());
      f1->SetParameters(5.e4,max_bin, 0.1);
      cout << endl << endl << "For Energy Range: " << Emin << " to " << Emax << endl;
      
      //sprintf(num, "his%i", i)     
      east->Fit("f1","RM");
      outfile0 << Emin << " " << Emax << " " << f1->GetParameter(1) << endl;
      delete f1;
      f1=NULL;

      TF1 *f2 = new TF1("f2", "landau", 0.05, 20.);
      //f2->SetParLimits(0, 0., 1.e10);
      //f2->SetParLimits(1, 0., 1.e7);

      max_bin = west->GetBinCenter(west->GetMaximumBin());
      f2->SetParameters(5.e4,max_bin,0.1);
      cout << endl << endl << "For Energy Range: " << Emin << " to " << Emax << endl;
      
      //sprintf(num, "his%i", i)     
      west->Fit("f2","RM");
      outfile1 << Emin << " " << Emax << " " << f2->GetParameter(1) << endl;
      

      //Energy_min=Energy_min + (double)binWidthEnergy;
      //Energy_max=Energy_max + (double)binWidthEnergy;
      
    }
  outfile0.close();
  outfile1.close();

}    

//TString inputsel = "&&nTaggableJetsPt20MatchedTagged>1";
void plotHhhTop(){
  gROOT->ProcessLine(".L tdrstyle.C");
  setTDRStyle();
  gROOT->SetBatch();

  //std::vector
  double binning_[]={0.,10.,20.,30.,40.,50.,60.,70.,80.,90.,100.,110.,120.,130.,140.,150.,160.,170.,180.,190.,200.,225.,250.,275.,300.,325.,350.,400.,500.,700.,1000.,1500.};

  TFile *HhhFile    = new TFile("/scratch/ojalvo/HhhGenStudies/Hhhttbb300.root");
  TFile *TOPFile = new TFile("/scratch/ojalvo/LTau_103/TOP2.root");
  TFile *MTEMFile = new TFile("/scratch/ojalvo/LTau_103/MTEM.root");

  TString selection="pt1>20&&pt2>20"+inputsel;//
  //TString selection="pt1>20&&pt2>20&&vertices>0&&diLeptons==0&&lPFIsoDB<0.1&&tauIso<1.5&&tightElectrons==0&&tightMuons<=1&&mt1<30&&J1Pt>20&&J2Pt>20";//

  TCanvas *Tcan= new TCanvas("Tcan","",100,20,600,800); Tcan->cd();  Tcan->SetFillColor(0); 

  TPad* pad1 = new TPad("pad1","The pad",0,0.2,1,1);
  TPad* pad2 = new TPad("pad2","the other pad",0,0,1,0.2);
  pad1->SetFillColor(0);
  //pad1->SetLogy();
  pad1->Draw();  pad1->cd();  pad1->SetLeftMargin(0.2);  pad1->SetBottomMargin(0.13);

  TLegend *leg = new TLegend(0.55,0.75,0.99,0.99);//  Leg = new TLegend(0.74,0.65,0.99,0.95);
  leg->SetFillColor(0);
  leg->SetFillStyle(0);
  leg->SetBorderSize(0);
  leg->SetHeader("");
  leg->SetShadowColor(0);

  TTree* HhhTree = (TTree*)HhhFile->Get("muTauEventTreeFinal/eventTree");

  //TH1F* Shape = new TH1F("hist1","hist1",31,&binning_[0]);
  TH1F* Shape = new TH1F("hist1","hist1",bins,lo,hi);
  HhhTree->Draw(var+">>+hist1","("+selection+")");
  Shape->SetLineColor(kBlue);  Shape->SetFillColor(kBlue);
  Shape->SetFillStyle(3002);
  Shape->SetLineWidth(4);
  //Shape->Sumw2();

  
  TTree* TOPTree = (TTree*)TOPFile->Get("muTauEventTreeFinal/eventTree");
  TH1F* Shape1 = new TH1F("hist2","hist2",bins,lo,hi);
  TOPTree->Draw(+var+">>+hist2","("+selection+")","same");

  Shape1->SetLineColor(kRed);
  Shape1->SetFillColor(kRed);
  Shape1->SetFillStyle(3002);
  Shape1->SetLineWidth(4);
  Shape1->Scale(Shape->Integral()/Shape1->Integral());

  //Shape1->Sumw2();

  TTree* MTEMTree = (TTree*)MTEMFile->Get("muTauEventTree/eventTree");
  TH1F* Shape2 = new TH1F("hist1","hist1",bins,lo,hi);
  MTEMTree->Draw(var+">>+hist1","("+selection+")","same");
  Shape2->SetLineColor(kBlack);
  Shape2->SetLineStyle(2);
  Shape2->SetLineWidth(4);
  Shape2->Scale(Shape->Integral()/Shape2->Integral());
  

  Shape->GetXaxis()->SetTitle(var);
  Shape->GetYaxis()->SetTitle("A.U.");
  Shape1->Draw("same");
  Shape->Draw("same");

  ///////////stupid thing to sent max bin
  float max, max1, max2 = 0;
  max = Shape->GetBinContent(Shape->GetMaximumBin());
  max1 = Shape1->GetBinContent(Shape1->GetMaximumBin());
  max2 = Shape2->GetBinContent(Shape2->GetMaximumBin());

  if(max>max1&&max>max2)
    Shape->SetMaximum(max*1.1);
  else if(max1>max&&max1>max2)
    Shape->SetMaximum(max1*1.1);
  else 
    Shape->SetMaximum(max2*1.1);

  ///
  leg->AddEntry(Shape,"300GeV H Signal","l");
  leg->AddEntry(Shape1,"t#bar{t}","l");
  leg->AddEntry(Shape2,"Z #rightarrow #tau#tau","l");
  //leg->AddEntry(Shape2,"Embedded 50 GeV<M(#mu#mu)<250 GeV","l");

  leg->Draw();

  
  Tcan->cd();
  /*
  pad2->Draw();
  pad2->cd();   
  pad2->SetBottomMargin(0.1);
  pad2->SetTopMargin(0.1);
  pad2->SetLeftMargin(0.2); 

  TH1F *Den = new TH1F("Den","hist",bins,lo,hi);
  TH1F *Num = new TH1F("Num","hist",bins,lo,hi);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      }
    }else{
      mchist[imc]->SetLineWidth(1);
      mchist[imc]->SetFillColor( colors[imc] );
      mchist[imc]->SetLineColor( 1 );
    }

    //mcstack->Add( mchist[imc] );

    if( !isSignal ){
      mcstack->Add( mchist[imc] );

      if( imc == 0 ){
	mctothist = (TH1F*) mchist[imc]->Clone();
	smtothist = (TH1F*) mchist[imc]->Clone();
      }
      else{
	mctothist->Add(mchist[imc]);
	smtothist->Add(mchist[imc]);
      }
    }
    else{
      mctothist->Add(mchist[imc]);
      sighist.push_back( mchist[imc] );
    }

    cout << "MC yield " << labels[imc] << " " << Form("%.2f",mchist[imc]->Integral()) << endl;
  }

  chdata->Draw(Form("TMath::Min(%s,%f)>>%s_datahist_%s",var,xmax-0.01,myvar,flavor),sel);

  if( overlayData ){

    float max = datahist->GetMaximum() + datahist->GetBinError(datahist->GetMaximumBin());
    if( mctothist->GetMaximum() > max ) max = mctothist->GetMaximum();
    if( log ) datahist->SetMaximum( 15 * max );
    else      datahist->SetMaximum( 1.4 * max );

    datahist->GetXaxis()->SetTitle(xtitle);
    datahist->Draw("E1");
    mcstack->Draw("samehist");
    
    for( unsigned int isig = 0 ; isig < sighist.size() ; isig++ ){
      sighist.at(isig)->Add(smtothist);
      sighist.at(isig)->Draw("samehist");
    } 

    datahist->Draw("sameE1");
    datahist->Draw("sameaxis");

    if(!log) datahist->GetYaxis()->SetRangeUser(0.,1.4*max);
    
  }
  else{
    float max = mctothist->GetMaximum();
    if( log ) mctothist->SetMaximum( 15 * max );
    else      mctothist->SetMaximum( 1.4 * max );

    mctothist->SetLineColor(0);
    mctothist->SetFillColor(0);

    mctothist->GetXaxis()->SetTitle(xtitle);
    mctothist->Draw("hist");
    mcstack->Draw("samehist");

    for( unsigned int isig = 0 ; isig < sighist.size() ; isig++ ){

void DrawQCDClosure(TString VAR,TString XTITLE)
{
  gROOT->ForceStyle();
  TString FileName[7] = {"QCD_HT200to300","QCD_HT300to500","QCD_HT500to700","QCD_HT700to1000","QCD_HT1000to1500","QCD_HT1500to2000","QCD_HT2000toInf"};
  float XSEC[7]       = {1.74e+6,3.67e+5,2.94e+4,6.524e+03,1.064e+03,121.5,2.542e+01};
  TFile *inf[7];
  TH1F  *h[7],*h1[7];

  TCanvas *can = new TCanvas("can_QCDClosure_"+VAR,"can_QCDClosure_"+VAR,900,600);
  can->cd(1);
  can->SetBottomMargin(0.3);
  can->SetRightMargin(0.15);

  for(int i=0;i<7;i++) {
    inf[i] = TFile::Open("Histo_"+FileName[i]+".root");
    TH1F *hTriggerPass = (TH1F*)inf[i]->Get("hadtopL/TriggerPass");
    h[i]   = (TH1F*)inf[i]->Get("hadtopL/h_"+VAR);
    h1[i]   = (TH1F*)inf[i]->Get("hadtop/h_"+VAR);
    h[i]->Sumw2();
    h1[i]->Sumw2();
    h[i]->Rebin(5);
    h1[i]->Rebin(5);
    h[i]->Scale(XSEC[i]/hTriggerPass->GetBinContent(1));
    h1[i]->Scale(XSEC[i]/hTriggerPass->GetBinContent(1));
    cout<<hTriggerPass->GetBinContent(1)<<endl;
  }

  TH1F *hQCD  = (TH1F*)h[0]->Clone("hQCD");
  TH1F *hQCD1 = (TH1F*)h1[0]->Clone("hQCD1");
  for(int i=0;i<7;i++) {
    hQCD->Add(h[i]);
    hQCD1->Add(h1[i]);
  } 
  hQCD->SetFillColor(kGray);

  hQCD->Scale(1./hQCD->Integral());
  hQCD1->Scale(1./hQCD1->Integral());
  hQCD->GetXaxis()->SetLabelSize(0.0);
  double max = 1.1*TMath::Max(hQCD->GetBinContent(hQCD->GetMaximumBin()),hQCD1->GetBinContent(hQCD1->GetMaximumBin()));
  hQCD->SetMinimum(1e-5);
  hQCD->SetMaximum(max);
  hQCD->Draw("hist");
  hQCD1->Draw("sameE");
  gPad->RedrawAxis();

  TLegend *leg = new TLegend(0.86,0.65,0.99,0.9);
  leg->SetFillColor(0);
  leg->SetTextFont(42);
  leg->SetTextSize(0.03);
  leg->AddEntry(hQCD,"Control","F");
  leg->AddEntry(hQCD1,"Signal","LP");
  leg->Draw();

  TH1F *hRatio = (TH1F*)hQCD1->Clone("Ratio");
  hRatio->Divide(hQCD);

  TPad* pad = new TPad("pad", "pad", 0., 0., 1., 1.);
  pad->SetTopMargin(0.7);
  pad->SetRightMargin(0.15);
  pad->SetFillColor(0);
  pad->SetFillStyle(0);
  pad->Draw();
  pad->cd(0);
  gPad->SetGridy();
  hRatio->GetXaxis()->SetTitle(XTITLE);
  hRatio->GetYaxis()->SetNdivisions(505);
  hRatio->GetYaxis()->SetRangeUser(0,2);
  hRatio->GetYaxis()->SetLabelSize(0.04);
  hRatio->Draw();
}

void DrawFits(TemplateFIT * FIT,FileSaver finalHistos,FileSaver Plots){

	std::string pathdata  = (FIT->GetName() + "/Fit Results/Data");
	std::string pathtemplP= (FIT->GetName() + "/Fit Results/ScaledTemplatesP");
	std::string pathtemplD= (FIT->GetName() + "/Fit Results/ScaledTemplatesD");	
	std::string pathtemplHe=(FIT->GetName() + "/Fit Results/ScaledTemplatesHe");	

	std::string pathfit   = (FIT->GetName() + "/Fit Results/FractionFits");
	std::string pathtrans = (FIT->GetName() + "/Fit Results/TrasnferFunctions");
	std::string pathres   = (FIT->GetName() + "/Fit Results/");
	
	TFile * infile = finalHistos.GetFile();	


	for(int i=1; i<FIT->GetBinning().size();i++){
		
		std::string pathbinP    = pathtemplP + "/Bin"+to_string(i);
		std::string pathbinD    = pathtemplD + "/Bin"+to_string(i);
		std::string pathbinHe   = pathtemplHe+ "/Bin"+to_string(i);
		std::string pathbindata = pathdata   + "/Bin"+to_string(i);
		std::string pathbinfit  = pathfit    + "/Bin"+to_string(i);

		std::vector<TH1F*> TemplatesP=GetListOfTemplates(infile, pathbinP);
		std::vector<TH1F*> TemplatesD=GetListOfTemplates(infile, pathbinD);		
		std::vector<TH1F*> TemplatesHe=GetListOfTemplates(infile, pathbinHe);		
		std::vector<TH1F*> Datas     =GetListOfTemplates(infile, pathbindata);
		std::vector<TH1F*> Fits      =GetListOfTemplates(infile, pathbinfit);
		std::vector<TH1F*> Transfer  =GetListOfTemplates(infile, pathtrans);

		cout<<pathbinHe<<" "<<TemplatesHe.size()<<" "<<TemplatesHe[0]<<endl;
	

		TCanvas * c1 = new TCanvas("Modified Templates");
		c1->SetCanvasSize(2000,1500);

		for(int j=TemplatesP.size()-1;j>=0;j--){
			if(j==0) PlotDistribution(gPad, TemplatesP[j],"Reconstructed Mass [GeV/c^2]","Counts",1,"same",1,TemplatesP[j]->GetBinContent(TemplatesP[j]->GetMaximumBin())*1.13,10);
			else     PlotDistribution(gPad, TemplatesP[j],"Reconstructed Mass [GeV/c^2]","Counts",colorbase + j,"same",1,TemplatesP[j]->GetBinContent(TemplatesP[j]->GetMaximumBin())*1.13,7,"",false,false,true);		
		}
	
		Plots.Add(c1);
		Plots.writeObjsInFolder((FIT->GetName()+"/Fits/Bin"+to_string(i)).c_str());	


		TCanvas * c2 = new TCanvas("Modified T. Fits");
                c2->SetCanvasSize(2000,1500);

		PlotDistribution(gPad, TemplatesP[0], "Reconstructed Mass [GeV/c^2]","Counts",2,"same",1,Datas[0]->GetBinContent(Datas[0]->GetMaximumBin())*1.13,10,"Original Protons MC Template");
		PlotDistribution(gPad, TemplatesD[0], "Reconstructed Mass [GeV/c^2]","Counts",4,"same",1,1e5,10,"Original Deuterons MC Template");
		PlotDistribution(gPad, TemplatesHe[0],"Reconstructed Mass [GeV/c^2]","Counts",3,"same",1,1e5,10,"Original He Fragm. MC Template");

		for(int j=TemplatesP.size()-1;j>=1;j--){
                        PlotDistribution(gPad, TemplatesP[j],"Reconstructed Mass [GeV/c^2]","Counts",2,"same",1,1e5,1,"",false,false,true);
			PlotDistribution(gPad, TemplatesD[j],"Reconstructed Mass [GeV/c^2]","Counts",4,"same",1,1e5,1,"",false,false,true);
                	//PlotDistribution(gPad, TemplatesHe[j],"Reconstructed Mass [GeV/c^2]","Counts",3,"same",1,1e5,1,"",false,false,true);
		}
		PlotDistribution(gPad, Datas[0],"Reconstructed Mass [GeV/c^2]","Counts",1,"ePsame",1,1e5,3,"ISS data",false,true);
	


		Plots.Add(c2);
                Plots.writeObjsInFolder((FIT->GetName()+"/Fits/Bin"+to_string(i)).c_str());
		
		
		TCanvas * c3 = new TCanvas("Template Fits");
                c3->SetCanvasSize(2000,1500);

		PlotDistribution(gPad, TemplatesP[0] ,"Reconstructed Mass [GeV/c^2]","Counts",2,"same",1,Datas[0]->GetBinContent(Datas[0]->GetMaximumBin())*1.13,2,"Original Protons MC Template");
		PlotDistribution(gPad, TemplatesD[0] ,"Reconstructed Mass [GeV/c^2]","Counts",4,"same",1,1e5,2,"Original Deuterons MC Template");
		PlotDistribution(gPad, TemplatesHe[0],"Reconstructed Mass [GeV/c^2]","Counts",3,"same",1,1e5,2,"Original He Fragm. MC Template");
		PlotDistribution(gPad, TemplatesP[1] ,"Reconstructed Mass [GeV/c^2]","Counts",2,"same",1,Datas[0]->GetBinContent(Datas[0]->GetMaximumBin())*1.13,10,"Best #chi^{2} Protons MC Template");
		PlotDistribution(gPad, TemplatesD[1] ,"Reconstructed Mass [GeV/c^2]","Counts",4,"same",1,1e5,10,"Best #chi^{2} Deuterons MC Template");
		PlotDistribution(gPad, TemplatesHe[1],"Reconstructed Mass [GeV/c^2]","Counts",3,"same",1,1e5,10,"Best #chi^{2} he Fragm. MC Template");
		

		PlotDistribution(gPad, Datas[0],"Reconstructed Mass [GeV/c^2]","Counts",1,"ePsame",1,1e5,3,"ISS data",false,true);
		if(Fits.size()>0) PlotDistribution(gPad, Fits[0],"Reconstructed Mass [gev/c^2]","counts",6,"same",1,1e5,4,"Fraction Fit");
	

		Plots.Add(c3);
                Plots.writeObjsInFolder((FIT->GetName()+"/Fits/Bin"+to_string(i)).c_str());

		TCanvas * c5 = new TCanvas("Transfer Functions");
                c5->SetCanvasSize(2000,1500);

		for(int j=0;j<Transfer.size();j++){
			Transfer[j]->Smooth(3);
			PlotDistribution(gPad, Transfer[j],"Reconstructed Mass [GeV/c^2]","Prim. / (Prim. + Sec.)",j,"hist,same",0,1,7,("Bin. "+to_string(j)).c_str(),false,false);
		}

		Plots.Add(c5);
                Plots.writeObjsInFolder((FIT->GetName()+"/Fits").c_str());

	

		TCanvas * c4 = new TCanvas("ChiSquare");
                c4->SetCanvasSize(2000,1500);
		gPad->SetLogz();	
		TH2F * Chi = (TH2F*) infile->Get((FIT->GetName()+"/Fit Results/Spreads/ChiSquare/ChiSquare Bin "+to_string(i)).c_str());
                Chi->GetZaxis()->SetRangeUser(0.2,100);
//.........这里部分代码省略.........