C++ TH2D::GetBinContent方法代码示例

void Occupancy::postProcessing()
{
  if (_postProcessed) return;

  // Generate the bounds of the 1D occupancy hist
  unsigned int totalHits = 0;
  unsigned int maxHits = 0;
  for (unsigned int nsens = 0; nsens < _device->getNumSensors(); nsens++)
  {
    Mechanics::Sensor* sensor = _device->getSensor(nsens);
    TH2D* occ = _hitOcc.at(nsens);

    for (unsigned int x = 0; x < sensor->getNumX(); x++)
    {
      for (unsigned int y = 0; y < sensor->getNumY(); y++)
      {
        const unsigned int numHits = occ->GetBinContent(x + 1, y + 1);
        totalHits += numHits;
        if (numHits > maxHits) maxHits = numHits;
      }
    }
  }

  TDirectory* plotDir = makeGetDirectory("Occupancy");

  std::stringstream name;
  std::stringstream title;

  name << "OccupancyDistribution";
  title << "Occupancy Distribution";

  _occDistribution = new TH1D(name.str().c_str(), title.str().c_str(),
                              100, 0, (double)maxHits / (double)totalHits);

  _occDistribution->SetDirectory(_dir);
  _occDistribution->GetXaxis()->SetTitle("Hits per trigger");
  _occDistribution->GetYaxis()->SetTitle("Pixels");
  _occDistribution->SetDirectory(plotDir);

  // Fill the occupancy distribution
  for (unsigned int nsens = 0; nsens < _device->getNumSensors(); nsens++)
  {
    Mechanics::Sensor* sensor = _device->getSensor(nsens);
    TH2D* occ = _hitOcc.at(nsens);

    for (unsigned int x = 0; x < sensor->getNumX(); x++)
    {
      for (unsigned int y = 0; y < sensor->getNumY(); y++)
      {
        const unsigned int numHits = occ->GetBinContent(x + 1, y + 1);
        _occDistribution->Fill((double)numHits / (double)totalHits);
      }
    }
  }

  _postProcessed = true;
}

TH2D* GetJetCorrFunc2D_ZYAM(int itrg, int jass)
{
                    TH2D* hcorr = (TH2D*)GetRawCorrFunc2D_ratio(itrg,jass);
                    TH2D* hcorr_clone = (TH2D*)hcorr->Clone(Form("corr_clone_itrg%d_jass%d",itrg,jass));
                    hcorr_clone->Reset();
		    for(int ietabin=1;ietabin<=hcorr->GetNbinsX();ietabin++)
		    {
                      TH1D* hcorrphi = (TH1D*)hcorr->ProjectionY(Form("corrphi_%d",ietabin),ietabin,ietabin,"e");  
                     // hcorrphi->SetAxisRange(0.35,1.57,"X");
                      double histminX = hcorrphi->GetBinCenter(hcorrphi->GetMinimumBin());
                      double histminY = hcorrphi->GetBinContent(hcorrphi->GetMinimumBin());
	   	      fitfunc = new TF1(Form("fitfunc_%d",ietabin),"[0]+[1]*(x-[2])*(x-[2])",0.6,1.2);		    
                      fitfunc->SetParameters(histminY,1,histminX);
                      fitfunc->SetParLimits(1,0,0.10000);
                  //    fitfunc->SetParLimits(2,0.35,1.57);
		      for(int ifit=0;ifit<3;ifit++) hcorrphi->Fit(Form("fitfunc_%d",ietabin),"RNO");                     		  
		      float level = fitfunc->GetParameter(0);
		      for(int iphibin=1;iphibin<=hcorr->GetNbinsY();iphibin++) 
                        hcorr_clone->SetBinContent(ietabin,iphibin,hcorr->GetBinContent(ietabin,iphibin)-level);
		    }
                    
                    float max = hcorr_clone->GetBinContent(hcorr_clone->GetMaximumBin());
                    hcorr_clone->SetAxisRange(ymin,max*1.3,"Z");		    
                    return hcorr_clone;
}

void runcorr(int filenum)
{


  string buffer;
  vector<string> listoffiles;
  int nlines = 0;
  ifstream infile("/net/hisrv0001/home/dav2105/corrana/makecorrhists/oldstuff/sortedforests.txt");

  if (!infile.is_open()) {
    cout << "Error opening file. Exiting." << endl;
    return;
  } else {
    while (!infile.eof()) {
      infile >> buffer;
      listoffiles.push_back(buffer);
      nlines++;
    }
  }


  corrana(listoffiles[filenum].data());
  
  double pttriglow[] = {1,4};
  double pttrighigh[] = {2,8};
  double ptasslow[] = {1,2};
  double ptasshigh[] = {2,3};
  int centmin[] = {0,4,8,12,16,20,24,28,32};
  int centmax[] = {40,8,12,16,20,24,28,32,36};
  TFile * outf = new TFile(Form("corrhists_%d.root",filenum),"recreate");
  
  // for(int i = 0 ; i < 9 ; ++i)
  // for(int i = 0 ; i < 2 ; ++i)
  for(int i = 0 ; i < 1 ; ++i)
  {
    cout<<"pt iteration "<<i<<endl;
    // for(int cent = 0 ; cent < 3 ; ++cent)
    // for(int cent = 0 ; cent < 9 ; ++cent)
    for(int cent = 0 ; cent < 1 ; ++cent)
    {
      cout<<"cent iteration "<<cent<<endl;
      TH2D * ttsig = TrackTrackSignal(pttriglow[i],pttrighigh[i],ptasslow[i],ptasshigh[i],centmin[cent],centmax[cent]);
      TH2D * ttbak = TrackTrackBackground(pttriglow[i],pttrighigh[i],ptasslow[i],ptasshigh[i],centmin[cent],centmax[cent]);
      // TCanvas * c1 = new TCanvas();
      // ttsig->Draw("surf1");
      // TCanvas * c2 = new TCanvas();
      // ttbak->Draw("surf1");
      // TCanvas * c3 = new TCanvas();
      TH2D * ttcorr = (TH2D*)ttsig->Clone(Form("corr_trg%d_%d_ass%d_%d_cmin%d_cmax%d",(int)pttriglow[i],(int)pttrighigh[i],(int)ptasslow[i],(int)ptasshigh[i],centmin[cent],centmax[cent]));
      ttcorr->Divide(ttbak);
      ttcorr->Scale(ttbak->GetBinContent(ttbak->FindBin(0,0)));
      ttcorr->GetXaxis()->SetRange(ttcorr->GetXaxis()->FindBin(-4.0),ttcorr->GetXaxis()->FindBin(4.0));
      ttcorr->GetYaxis()->SetRange(ttcorr->GetYaxis()->FindBin(-3.1415926/2.0),ttcorr->GetYaxis()->FindBin(3*3.1415926/2.0));
      // ttcorr->Draw("surf1");
    }
  }
  
  outf->Write();
  outf->Close();
}

void fukuchiPhantom2microPhantom_muon2(){

  gROOT->Reset();
  TFile *fin = new TFile("fukuchiPhantom.root");
  TH2D *hin = (TH2D *)fin->Get("phantom");

  int nbinx = hin->GetXaxis()->GetNbins();
  int nbiny = hin->GetYaxis()->GetNbins();


  TH2D *hout = new TH2D("phantom","phantom",nbinx,-0.75,0.75,nbiny,-0.75,0.75);

  for(int i=1;i<=nbinx;i++){
    for(int j=1;j<=nbiny;j++){
      double v = hin->GetBinContent(i,j);
      hout->SetBinContent(i,j,v);
    }
  }


  TFile *fout = new TFile("microPhantom_muon2.root","recreate");
  fout->cd();
  hout->Write();
  fout->Close();

}

void HitInfo::postProcessing()
{
  for (unsigned int nsens = 0; nsens < _device->getNumSensors(); nsens++)
  {
    TH2D* map = _totMap.at(nsens);
    TH2D* count = _totMapCnt.at(nsens);

    for (Int_t x = 1; x <= map->GetNbinsX(); x++)
    {
      for (Int_t y = 1; y <= map->GetNbinsY(); y++)
      {
        const double average = map->GetBinContent(x, y) / count->GetBinContent(x, y);
        map->SetBinContent(x, y, average);
      }
    }
  }
}

void Matching::postProcessing()
{
  for (unsigned int nsens = 0; nsens < _dutDevice->getNumSensors(); nsens++)
  {
    Mechanics::Sensor* sensor = _dutDevice->getSensor(nsens);

    TH2D* tot = _inPixelTot.at(nsens);
    TH2D* hits = _inPixelTracks.at(nsens);

    assert(hits->GetNbinsX() == tot->GetNbinsX() &&
           hits->GetNbinsY() == tot->GetNbinsY() &&
           "Matching: histograms should have the same number of bins");

    for (int nx = 1; nx <= tot->GetNbinsX(); nx++)
    {
      for (int ny = 1; ny <= tot->GetNbinsY(); ny++)
      {
        const unsigned int num = hits->GetBinContent(nx, ny);
        if (num == 0) continue;
        const double average = tot->GetBinContent(nx, ny) / (double)num;
        tot->SetBinContent(nx, ny, average);
      }
    }

    for (unsigned int axis = 0; axis < 2; axis++)
    {
      TH1D* dist = axis ? _matchDistX.at(nsens) :
                          _matchDistY.at(nsens);

      double pixelWidth = axis ? sensor->getPitchX() :
                                 sensor->getPitchY();
      double beamSigma = axis ? _refDevice->getSensor(0)->getPitchX() :
                                _refDevice->getSensor(0)->getPitchY();

      TF1* fit = Processors::fitPixelBeam(dist, pixelWidth, beamSigma, false);

      std::stringstream ss;
      ss << _dutDevice->getName() << sensor->getName() << "PixelBeamFit"
         << ( axis ? "X" : "Y") << _nameSuffix;
      _plotDir->WriteObject(fit, ss.str().c_str());
    }
  }
}

void DrawTH2DZ(TFile* File, char* Histos_Name, char* legend, char* Title, char* Xlegend, char* Ylegend, double xmin, double xmax, double ymin, double ymax,  bool save, char* save_path, bool Normalization)
{
   TH2D* Histos = File->Get(Histos_Name);
   Histos->SetTitle();
   Histos->SetStats(kFALSE);
   Histos->GetXaxis()->SetTitle(Xlegend);
   Histos->GetYaxis()->SetTitle(Ylegend);
   Histos->GetYaxis()->SetTitleOffset(1.20);
//   Histos->GetYaxis()->SetTitleOffset(1.0);
   if(xmin!=xmax)Histos->SetAxisRange(xmin,xmax,"X");
   if(ymin!=ymax)Histos->SetAxisRange(ymin,ymax,"Y");

   if(Normalization){
      for(int x=0;x<Histos->GetXaxis()->GetNbins();x++){     
         TH1D* tmp = Histos->ProjectionY("",x,x);
	 double Integral = tmp->Integral();
         if(Integral==0)continue;
         double Factor = 1/Integral;
         for(int y=0;y<Histos->GetYaxis()->GetNbins();y++){
            Histos->SetBinContent(x,y, Histos->GetBinContent(x,y)*Factor );
            Histos->SetBinError  (x,y, Histos->GetBinError  (x,y)*Factor );
         }
      }
   }




   Histos->Draw("COLZ");
   gPad->SetLogz(1);
/*   c1->Update();
   TPaletteAxis* palette = (TPaletteAxis*)Histos->GetListOfFunctions()->FindObject("palette");
   palette->SetLabelOffset(0.1);
   palette->SetTitleOffset(0.1);
   c1->Modified();
*/

  
   if(save==1){
//        char path[255]; sprintf(path,"Pictures/PNG/%s.png",save_path);  c1->SaveAs(path);
//        char path[255]; sprintf(path,"Pictures/EPS/%s.eps",save_path);  c1->SaveAs(path);
//        char path[255]; sprintf(path,"Pictures/C/%s.C"  ,save_path);  c1->SaveAs(path);
           c1->SaveAs(save_path);

   }
}

void MakeGraph()
{
    InitPar();

    for(Int_t i = 0; i < kNData; i++) {
        TFile *file = new TFile(Form("Att%d.root", i * 100), "OPEN");
        TH2D *hist = (TH2D*)file->Get("HisPlate");

        thickness[i] = i * 100.;
        Double_t nHit[16] = {0};

        for(Int_t x = 0; x < 12; x++) {
            for(Int_t y = 0; y < 8; y++) {
                Int_t index = GetIndex(GetDistance(x - 5.5, y - 3.5));
                if(index == -1) cout << "Error" << endl;
                else {
                    nHit[index]++;
                    depEnergy[index][i] += hist->GetBinContent(x + 1, y + 1);
                }
            }
        }

        for(Int_t index = 0; index < 16; index++) {
            if( nHit[index > 0]) depEnergy[index][i] /= nHit[index];
            else cout << "Really?" << endl;
        }

        delete hist;
        delete file;
    }

    for(Int_t index = 0; index < 16; index++) {
        Double_t max =  depEnergy[index][0];
        for(Int_t iT = 0; iT < 31; iT++) {
            depEnergy[index][iT] /= max;
        }
    }

    for(Int_t i = 0; i < 16; i++) {
        gr[i] = new TGraph(kNData, thickness, depEnergy[i]);
        fitFnc[i] = new TF1(Form("fitFnc%02d", i), "pol9", 0, 3000);
        fitFnc[i]->FixParameter(0, 1.);
        gr[i]->Fit(fitFnc[i], "R");
    }

}

void mytestreweight(){


  TFile*finput=new TFile("fworking.root");
  TH1D*hBPtFONLL=(TH1D*)finput->Get("hBPt");
  TH1D*hDPtFONLL=(TH1D*)finput->Get("hDPt");
  TH2D*hD=(TH2D*)finput->Get("hD");
  
   hD->SetXTitle("B p_{T} (GeV/c)");
   hD->SetYTitle("D^{0} p_{T} (GeV/c)");
   hBPtFONLL->SetXTitle("B p_{T} (GeV/c)");
   hDPtFONLL->SetXTitle("D p_{T} (GeV/c)");
  
  TH2D *hDreweight = (TH2D*)hD->Clone("hDreweight");
  TH1D *hBPtPythia = hD->ProjectionX("hBPtPythia");
  TH1D *hDPtPythia = hD->ProjectionY("hDPtPythia");

   for (int x=1;x<=hDreweight->GetNbinsX()+1;x++){  //loop over the B pt bins
      if (hBPtPythia->GetBinContent(x)==0) continue;
      double ratio = hBPtFONLL->GetBinContent(x)/hBPtPythia->GetBinContent(x);  // in each pt bin we calculate the ratio of pythiaB/FONLLB 
       
      for (int y=1;y<=hDreweight->GetNbinsY()+1;y++){   //loop over the D pt bins
         double ratio2 = ratio;                                              //copy the ratio for each B bin
         double val = hDreweight->GetBinContent(x,y)*ratio2;
         double valError = hDreweight->GetBinError(x,y)*ratio2;
	 hDreweight->SetBinContent(x,y,val);
	 hDreweight->SetBinError(x,y,valError);
      }   
   }

  TH1D *hBPtPythiaReweight = hDreweight->ProjectionX("hBPtPythiaReweight");
  TH1D *hDPtPythiaReweight = hDreweight->ProjectionY("hDPtPythiaReweight");
  
  hBPtPythiaReweight->Divide(hBPtFONLL);
   
  TCanvas*canvas=new TCanvas("canvas","canvas",1000,500);  
  canvas->SetLogy();
  canvas->Divide(3,1);
  canvas->cd(1);
  hBPtFONLL->Draw();
  canvas->cd(2);
  hBPtPythia->Draw();
  canvas->cd(3);
  hBPtPythiaReweight->Draw();

}

std::pair<Matrix, Matrix> roothist_to_matrix(const TH2D & hist, bool transpose){
    int ngen = hist.GetNbinsY();
    int nreco = hist.GetNbinsX();
    if(transpose){
        swap(ngen, nreco);
    }
    
    Matrix m = Matrix(nreco, ngen);
    Matrix m_e = Matrix(nreco, ngen);
    
    for(int i=0; i<nreco; ++i){
        for(int j=0; j<ngen; ++j){
            int nx = i+1;
            int ny = j+1;
            if(transpose) swap(nx, ny);
            m(i,j) = hist.GetBinContent(nx, ny);
            m_e(i,j) = hist.GetBinError(nx, ny);
        }
    }
    return pair<Matrix, Matrix>(move(m), move(m_e));
}

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

  gDirectory->ReadAll(); // load histos

  TList * lst = gDirectory->GetList();
  cout << lst->GetName() << endl;
  cout << lst->GetTitle() << endl;
  cout << "size    " << lst->GetSize() << endl;
  cout << "entries " << lst->GetEntries() << endl;
  cout << "last    " << lst->LastIndex() << endl;

  TIterator *iter = lst->MakeIterator();
  int ii = 0;
  TObject *obj;
  TH1D *h;
  TH1D *h0;
  TH2D *H;
  TH2D *H0;

  while( obj = iter->Next() ){
    ii++;
    cout << setw(4) << ii << ": ";
    cout << obj->ClassName() << " ";
    cout << obj->InheritsFrom("TH1D") << " ";
    cout << obj->GetName() << " \"";
    cout << obj->GetTitle() << "\"";
    cout << endl;
    //    if( obj->ClassName() == "TH1D" ){
    if( obj->InheritsFrom("TH1D") ){
      h = (TH1D*) obj;
      cout << "       1D";
      cout << h->GetNbinsX() << " bins, ";
      cout << h->GetEntries() << " entries, ";
      cout << h->GetSumOfWeights() << " inside, ";
      cout << h->GetBinContent(0) << " under, ";
      cout << h->GetBinContent(h->GetNbinsX()+1) << " over";
      cout << endl;

      f0.cd(); // output file

      //      TH1D* h0 = (TH1D*) h->Clone();
      h0 = h; // copy
      h0->Write(); // write to file f0
      
      f1.cd(); // back to file 1 for the loop
    }
    else{

      if( obj->InheritsFrom("TH2D") ){

	H = (TH2D*) obj;
	cout << "       2D";
	cout << H->GetNbinsX() << " bins, ";
	cout << H->GetEntries() << " entries, ";
	cout << H->GetSumOfWeights() << " inside, ";
	cout << H->GetBinContent(0) << " under, ";
	cout << H->GetBinContent(H->GetNbinsX()+1) << " over";
	cout << endl;
	
	f0.cd(); // output file
	
	H0 = H; // copy
	H0->Write(); // write to file f0
	
	f1.cd(); // back to file 1 for the loop
      }
      else cout << "other class " << obj->ClassName() << endl;

/** 
 * 
 * 
 * @param o 
 * @param useWeights 
 * @param correct
 *
 * @ingroup pwglf_forward_scripts_tests
 */
void
TestPoisson(Double_t o=.3, bool useWeights=false, bool correct=true)
{
  const char* load = "$ALICE_PHYSICS/PWGLF/FORWARD/analysis2/scripts/LoadLibs.C";
  if (!gROOT->GetClass("AliAODForwardMult")) {
    gROOT->Macro(load);
    gROOT->GetInterpreter()->UnloadFile(gSystem->ExpandPathName(load));
  }
  
  // --- Parameters of this script -----------------------------------
  Int_t      nBin =  5;  // Our detector matrix size 
  Int_t      nMax = TMath::Max(Int_t(nBin * nBin * o + .5)+nBin/2,nBin);  
  Int_t      nEv  = 10000; // Number of events
  Double_t   mp   = o;   // The 'hit' probability 


  TH2D* base = new TH2D("base", "Basic histogram", 
			nBin,-.5, nBin-.5, nBin, -.5, nBin-.5);
  base->SetXTitle("#eta");
  base->SetYTitle("#varphi");
  base->SetDirectory(0);
  base->SetOption("colz");

  Int_t tN1=nMax;    Double_t tMin1; Double_t tMax1;
  Int_t tN2=nMax*10; Double_t tMin2; Double_t tMax2=nMax;
  MakeIntegerAxis(tN1, tMin1, tMax1);
  MakeIntegerAxis(tN2, tMin2, tMax2);
  TH2D* corr = new TH2D("comp", "Comparison", 
			tN1, tMin1, tMax1, tN2, tMin2, tMax2);
  corr->SetXTitle("Input");
  corr->SetYTitle("Poisson");
  corr->SetDirectory(0);
  corr->SetOption("colz");
  corr->SetStats(0);
  TLine* lcorr = new TLine(0, 0, tMax2, tMax2);

  Int_t mm = TMath::Max(Int_t(nBin * o + .5),nBin/2);
  tN2=mm*10; tMax2 = mm;
  MakeIntegerAxis(tN2, tMin2, tMax2);
  Info("", "Making mean w/nbins=%d,range=[%f,%f]", tN2, tMin2, tMax2);
  TH2D* mean = new TH2D("mean", "Mean comparison", 
			tN2, tMin2, tMax2, tN2, tMin2, tMax2);
  mean->SetXTitle("Input");
  mean->SetYTitle("Poisson");
  mean->SetDirectory(0);
  mean->SetOption("colz");
  mean->SetStats(0);
  TLine* lmean = new TLine(tMin2, tMin2, tMax2, tMax2);

  TH1D* dist = new TH1D("dist", "Distribution of hits", tN1, tMin1, tMax1);
  dist->SetXTitle("s");
  dist->SetYTitle("P(s)");
  dist->SetFillColor(kRed+1);
  dist->SetFillStyle(3001);
  dist->SetDirectory(0);

  TH1D* diff = new TH1D("diff", "P-T", 100, -25, 25);
  diff->SetXTitle("Difference");
  diff->SetFillColor(kRed+1);
  diff->SetFillStyle(3001);
  diff->SetYTitle("Prob");

  AliPoissonCalculator* c = new AliPoissonCalculator("ignored");
  c->Init(nBin ,nBin);

  for (Int_t i = 0; i < nEv; i++) { 
    c->Reset(base);
    base->Reset();

    for (Int_t iEta = 0; iEta < nBin; iEta++) { 
      for (Int_t iPhi = 0; iPhi < nBin; iPhi++) { 
	// Throw a die 
	Int_t m = gRandom->Poisson(mp);
	dist->Fill(m);

	// Fill into our base histogram 
	base->Fill(iEta, iPhi, m);

	// Fill into poisson calculator 
	c->Fill(iEta, iPhi, m > 0, (useWeights ? m : 1));
      }
    }
    // Calculate the result 
    TH2D* res = c->Result(correct);
    
    // Now loop and compare 
    Double_t mBase = 0;
    Double_t mPois = 0;
    for (Int_t iEta = 0; iEta < nBin; iEta++) { 
      for (Int_t iPhi = 0; iPhi < nBin; iPhi++) { 
	Double_t p = res->GetBinContent(iEta, iPhi);
//.........这里部分代码省略.........

void Plot_2D(string filename){
  
  //string filename = "data/scan_CH1-50_masked.root";
 
  // string varXname = "VMM # [1-8]";
  // string varYname = "CH # [1-64]";

  string varXname = "Pulsed CH # [1-64]";
  string varYname = "Recorded CH # [1-64]";
  string varZname = "# Recorded Pulses [/100 pulsed]";
  //string varZname = "#sigma(PDO) / #bar{PDO}";
  // int Nx = 8;
  // double Xmin = 0.5;
  // double Xmax = 8.5;
  int Nx = 64;
  double Xmin = 0.5;
  double Xmax = 64.5;
  // int Ny = 64;
  // double Ymin = 0.5;
  // double Ymax = 64.5;
  int Ny = 4097;
  double Ymin = -0.5;
  double Ymax = 4096.5;

  ///////////////////////////////////////////////////////
  setstyle(0);
  
  TChain* tree = new TChain("VMM_data","VMM_data");

  tree->AddFile(filename.c_str());

  MMFE8Base* base = new MMFE8Base(tree);

  int N = tree->GetEntries();

  TH2D* hist = new TH2D("hist","hist",
			Nx, Xmin, Xmax,
			Ny, Ymin, Ymax);

  TH2D* hist2 = new TH2D("hist2","hist2",
			 Nx, Xmin, Xmax,
			 Ny, Ymin, Ymax);

  TH2D* histN = new TH2D("histN","histN",
			 Nx, Xmin, Xmax,
			 Ny, Ymin, Ymax);

  for(int i = 0; i < N; i++){
    base->GetEntry(i);

    if(base->VMM != 0)
      continue;

    if(base->CHword == 24)
      continue;

    if(base->Delay != 30)
      continue;
    
    // hist->Fill(base->CHpulse,base->CHword,base->PDO);
    //histN->Fill(base->CHpulse,base->CHword);
    histN->Fill(base->CHpulse,base->BCID);

    // if(base->CHpulse != 5)
    //   continue;
    // base->CHpulse = base->CHword;

    hist->Fill(base->VMM,base->CHpulse,base->PDO);
    hist2->Fill(base->VMM,base->CHpulse,base->PDO*base->PDO);
    // if(histN->GetBinContent(base->VMM,base->CHpulse) <= 0)
    //   histN->SetBinContent(base->VMM,base->CHpulse,base->BCID);
    // else
    //   if(histN->GetBinContent(base->VMM,base->CHpulse) > base->BCID)
    // histN->SetBinContent(base->VMM,base->CHpulse,base->BCID);
    // histN->Fill(base->VMM,base->CHpulse);
  
  }

  for(int x = 0; x < Nx; x++){
    for(int y = 0; y < Ny; y++){
      double v = hist->GetBinContent(x+1,y+1);
      double v2 = hist2->GetBinContent(x+1,y+1);
      double N = histN->GetBinContent(x+1,y+1);
      double vbar = v/max(int(N),1);
      double v2bar = v2/max(int(N),1);
      //hist->SetBinContent(x+1,y+1,sqrt(v2bar-vbar*vbar)/vbar);
      hist->SetBinContent(x+1,y+1,vbar);
    }
  }

  TCanvas* can = new TCanvas("can","can",600,500);
  can->SetLeftMargin(0.15);
  can->SetRightMargin(0.22);
  can->SetBottomMargin(0.15);
  can->SetTopMargin(0.08);

  can->Draw();
  can->SetGridx();
  can->SetGridy();
  
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    canvas->Divide( kNChan, 2 );

    // ok now start the loop on channels
    for ( size_t iChan = 0 ; iChan < kNChan ; ++iChan ) { 
	
      if ( iDetector == 0 ) channel[iChan] = iChan - 0.5;

      string tempName = "NoiseMap_d" + toString( iDetector ) + "_l" + toString( loop )	+ "_ch" + toString( iChan ) ;
      string tempTitle = "NoiseMap Det. " + toString( iDetector ) + " - Ch. " + toString( iChan ) ;

      TH2D * noiseMapCh = new TH2D ( tempName.c_str() , tempTitle.c_str(), 
				     kXPixel / kNChan , -0.5 + xLimit[ iChan ] , -0.5 + xLimit[ iChan + 1 ],
				     kYPixel, -0.5, -0.5 + kYPixel );
      noiseMapCh->SetXTitle("X [pixel]");
      noiseMapCh->SetYTitle("Y [pixel]");
      noiseMapCh->SetZTitle("Noise [ADC]");
      noiseMapCh->SetStats( false );
      outputHistoList->Add( noiseMapCh ) ;

      tempName = "NoiseDist_d" + toString( iDetector ) + "_l" + toString( loop )	+ "_ch" + toString( iChan ) ;
      
      tempTitle = "NoiseDist Det. " + toString( iDetector ) + " - Ch. " + toString( iChan ) ; 

      TH1D * noiseDistCh = new TH1D( tempName.c_str(), tempTitle.c_str(), 50, 0., 10. );
      noiseDistCh->SetXTitle("Noise [ADC]");
      noiseDistCh->SetLineColor( kColor[iDetector]  );
      noiseDistCh->SetLineStyle( iChan + 2 );
      noiseDistCh->SetLineWidth( 2 );
      outputHistoList->Add( noiseDistCh );

      // let's start looping on pixels now
      for ( size_t yPixel = 1 ; yPixel <= kYPixel ; ++yPixel ) {
	for ( size_t xPixel = xLimit[ iChan ] + 1; xPixel <= xLimit[ iChan +1 ] ; ++xPixel ) {
	  double noise = noiseMap->GetBinContent( xPixel , yPixel );
	  noiseMapCh->Fill( xPixel - 1 , yPixel - 1, noise );
	  noiseDistCh->Fill( noise );
	  
	}
      }

      canvas->cd( iChan + 1 ) ;
      noiseMapCh->Draw("colz");
      canvas->cd( iChan + kNChan + 1  );
      noiseDistCh->Draw();
      
      topPad->cd( iDetector + 1 );
      if ( iChan == 0 ) {
	noiseDistCh->Draw();
      } else {
	noiseDistCh->Draw("same");
      }

      middlePad->cd( iChan + 1 );
      if ( iDetector == 0 ) {
	noiseDistCh->Draw();
      } else {
	noiseDistCh->Draw("same");
      }


      noiseMean[ kNChan * iDetector + iChan ] = noiseDistCh->GetMean();
      noiseRMS[ kNChan * iDetector  + iChan ] = noiseDistCh->GetRMS();
    }
    canvas->Write();

  }

//.........这里部分代码省略.........
    hAllFound->Draw("same");  // MC pt distribution
  }
 
 
  /*

  .L ~/ITSupgrade/BuildDetector/DetectorK.cxx+
  
  // All NEW
  DetectorK its("ALICE","ITS");
  its.MakeAliceAllNew(0);
  its.SetMaxRadiusOfSlowDetectors(0.01);
  its.SolveViaBilloir(0);
  TGraph *c = its.GetGraphRecoEfficiency(0,3,2);
  c->Draw("C");


  // Current
  DetectorK its("ALICE","ITS");
  its.MakeAliceCurrent(0,0);
  its.SetMaxRadiusOfSlowDetectors(0.01);
  its.SolveViaBilloir(0);
  TGraph *c = its.GetGraphRecoEfficiency(0,4,2);
  c->Draw("C");

  */

  TCanvas *c3 = new TCanvas("c3","impact");//,200,10,900,900);
  c3->Divide(2,1); c3->cd(1);
  // Impact parameter

  // Impact parameter resolution ---------------
  h2Ddca->Draw("colz");
  h2Ddca->FitSlicesY() ;
  TH2D *dcaM = (TH2D*)gDirectory->Get("dca2D_1"); dcaM->Draw("same");
  TH2D *dcaRMS = (TH2D*)gDirectory->Get("dca2D_2"); //dcaRMS->Draw();
  TGraphErrors *d0 = new TGraphErrors(); 
  for (Int_t ibin =1; ibin<=dcaRMS->GetXaxis()->GetNbins(); ibin++) {
    d0->SetPoint(     ibin-1,dcaRMS->GetBinCenter(ibin),dcaRMS->GetBinContent(ibin)*1e4); // microns
    d0->SetPointError(ibin-1,0,dcaRMS->GetBinError(ibin)*1e4); // microns
  }
  d0->SetMarkerStyle(21);
  d0->SetMaximum(200);  d0->SetMinimum(0);
  d0->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)");
  d0->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)");
  d0->SetName("dca");  d0->SetTitle("DCAvsPt");

  c3->cd(1);  h2Ddca->Draw("surf2");
  c3->cd(2);  d0->Draw("APE");

  // PT RESOLUTION ------------
  TCanvas *c4 = new TCanvas("c4","pt resolution");//,200,10,900,900);
  c4->Divide(2,1); c4->cd(1);
  // Impact parameter
  h2Dpt->Draw("colz");
  h2Dpt->FitSlicesY() ;
  TH2D *dPtM = (TH2D*)gDirectory->Get("dPt2D_1"); dPtM->Draw("same");
  TH2D *dPtRMS = (TH2D*)gDirectory->Get("dPt2D_2"); // dPtRMS->Draw("");
  TGraphErrors *gPt = new TGraphErrors(); 
  for (Int_t ibin =1; ibin<=dPtRMS->GetXaxis()->GetNbins(); ibin++) {
    gPt->SetPoint(     ibin-1,dPtRMS->GetBinCenter(ibin),dPtRMS->GetBinContent(ibin)); 
    gPt->SetPointError(ibin-1,0,dPtRMS->GetBinError(ibin)); 
  }
  gPt->SetMarkerStyle(21);
  gPt->SetMaximum(20);  gPt->SetMinimum(0);
  gPt->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)");
  gPt->GetYaxis()->SetTitle("relative momentum resolution (%)");
  gPt->SetName("dPt");  gPt->SetTitle("DPTvsPt");

  c4->cd(1);  h2Dpt->Draw("surf2");
  c4->cd(2);  gPt->Draw("APE");


  // EXPORT --------

  TFile f("histos.root","RECREATE");

  hMultCount->Write();
  hAllMC->Write();
  hAllFound->Write();
  hImperfect->Write();
  hPerfect->Write();
  hNoMCTrack->Write();

  hPurity->Write();
  hEff->Write();
  hFake->Write();
  hAnna->Write();

  h2Ddca->Write();
  d0->Write();

  h2Dpt->Write();
  gPt->Write();

  f.Close();

  return;

}