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

void fadc_fit_heights(std::string fadc_id,std::string hist_type)
{
  /*****************************************************************/
  // Prepare the canvas
  gStyle->SetOptFit(1111);
  TCanvas *AlCapCanvas = (TCanvas *) gROOT->GetListOfCanvases()->At(0);
  AlCapCanvas->Clear();
  AlCapCanvas->Divide(4,2);

  //  gROOT->ProcessLine(".L common/get_histogram.C"); // get_histogram() is called here
  /*****************************************************************/
  const int n_channels = 8;
  std::string bank_names[n_channels] = {"Na", "Nb", "Nc", "Nd", "Ne", "Nf", "Ng", "Nh"};
  std::string name;

  double mean = 0;
  double rms = 0;

  for (int iChn = 0; iChn < n_channels; iChn++) {
    name=bank_names[iChn]+fadc_id;
    TH1* hist = get_histogram(name, hist_type);
    mean = hist->GetMean();
    rms = hist->GetRMS();

      AlCapCanvas->cd(iChn+1);

    if (hist->Integral()!=0) {
      TF1* gaus = new TF1("gaus","gaus",mean-rms,mean+rms);
      hist->Fit("gaus","Q");
      hist->GetXaxis()->SetRange(mean-2*rms,mean+2*rms);
      hist->Draw();
    }
  }
}

//________________________________________________________
void GFOverlay::CreateFillMeanRms(const TObjArray &hists, Int_t layer, const char *dirName,
				  std::vector<TH1*> &meanHists, std::vector<TH1*> &rmsHists) const
{
  // fill mean/rms from hists into the corresponding meanHists/rmsHists
  // if these are empty, create one hist for each slot of hists (even for empty ones!)
  if (hists.IsEmpty()) return;
  TH1 *h1 = 0;
  for (Int_t iH = 0; !h1 && iH < hists.GetEntriesFast(); ++iH) {
    h1 = static_cast<TH1*>(hists[iH]);
  }
  if (!h1 || h1->GetDimension() > 1) return; // only for 1D hists
  
  if (meanHists.empty()) { // create mean/RMS hists if not yet done
    const Float_t min = h1->GetXaxis()->GetXmin()/3.;
    const Float_t max = h1->GetXaxis()->GetXmax()/3.;
    const Int_t nBins = h1->GetNbinsX()/2;
    for (Int_t iHist = 0; iHist < hists.GetEntriesFast(); ++iHist) {
      TH1 *hMean = new TH1F(Form("mean%d_%d", layer, iHist), Form("%s: mean", dirName),
			    nBins, min, max);
      meanHists.push_back(hMean);
      TH1 *hRms = new TH1F(Form("rms%d_%d", layer, iHist), Form("%s: RMS", dirName),
			   nBins, 0., max);
      rmsHists.push_back(hRms);
    }
  }

  // now fill mean and rms hists
  for (Int_t iHist = 0; iHist < hists.GetEntriesFast(); ++iHist) {
    TH1 *h = static_cast<TH1*>(hists[iHist]);
    if (!h) continue;
    meanHists[iHist]->Fill(h->GetMean());
    rmsHists[iHist]->Fill(h->GetRMS());
  }
}

  void getPlotData() {
    TH1 * h = (TH1*) m_file->Get(m_direc.c_str());

    for (int i=0; i<h->GetXaxis()->GetNbins(); i++) {
      m_xs.push_back(h->GetXaxis()->GetBinCenter(i));
      m_ys.push_back(h->GetBinContent(i));
    }

    m_plot->m_xAxisTitle = std::string(h->GetXaxis()->GetTitle());
    m_plot->m_yAxisTitle = std::string(h->GetYaxis()->GetTitle());
    m_plot->m_title = std::string(h->GetTitle());

    std::stringstream ssN, ssMu, ssSig, ssUF, ssOF;
    ssN << std::setprecision(4) << h->GetEntries();
    ssMu << std::setprecision(4) << h->GetMean();
    ssSig << std::setprecision(4) << h->GetRMS();
    ssUF << std::setprecision(4) << h->GetBinContent(0);
    ssOF << std::setprecision(4) << h->GetBinContent(h->GetNbinsX() + 1);

    m_statsTitles.push_back("N:");
    m_statsTitles.push_back("mu:");
    m_statsTitles.push_back("sig:");
    m_statsTitles.push_back("UF:");
    m_statsTitles.push_back("OF:");

    m_statsValues.push_back(ssN.str());
    m_statsValues.push_back(ssMu.str());
    m_statsValues.push_back(ssSig.str());
    m_statsValues.push_back(ssUF.str());
    m_statsValues.push_back(ssOF.str());
  }

void DrawStatBox(TObject** Histos, std::vector<char*> legend, bool Mean, double X, double Y, double W, double H)
{
   int    N             = legend.size();
   char   buffer[255];

   if(Mean)H*=3;
   for(int i=0;i<N;i++){
           TPaveText* stat = new TPaveText(X,Y-(i*H), X+W, Y-(i+1)*H, "NDC");
	   TH1* Histo = (TH1*)Histos[i];
           sprintf(buffer,"Entries : %i\n",(int)Histo->GetEntries());
           stat->AddText(buffer);

           if(Mean){
           sprintf(buffer,"Mean    : %6.2f\n",Histo->GetMean());
           stat->AddText(buffer);

           sprintf(buffer,"RMS     : %6.2f\n",Histo->GetRMS());
           stat->AddText(buffer);
           }

           stat->SetFillColor(0);
           stat->SetLineColor(Color[i]);
           stat->SetTextColor(Color[i]);
           stat->SetBorderSize(0);
           stat->SetMargin(0.05);
           stat->SetTextAlign(12);
           stat->Draw();
   }
}

void logStatisticsPar(std::ostream& out, RooDataSet *dataSet, RooRealVar *realVar, int nBins, double chi2, const RooArgList &variables)
{
    TH1 *histogram = dataSet->createHistogram(Form("h%s", dataSet->GetName()), *realVar, RooFit::Binning(nBins));
    
    // Create the TeX file
    out << "\\documentclass[10pt]{article}" << std::endl;
    out << "\\usepackage[usenames]{color} %used for font color" << std::endl;
    out << "\\usepackage{fontspec}" << std::endl;
    out << "\\usepackage{xunicode}" << std::endl;
    out << "\\usepackage{xltxtra}" << std::endl;
    out << "\\defaultfontfeatures{Scale=MatchLowercase}" << std::endl;
    out << "\\setromanfont[Mapping=tex-text]{Myriad Pro}" << std::endl;
    out << "\\setsansfont[Mapping=tex-text]{Myriad Pro}" << std::endl;
    out << "\\setmonofont{Monaco}" << std::endl;
    out << "\\begin{document}" << std::endl;
    out << "\\thispagestyle{empty}" << std::endl;
    out << "\\setlength{\\tabcolsep}{1ex}" << std::endl;
    out << "\\setlength{\\fboxsep}{0ex}" << std::endl;
    out << "{\\fontsize{7pt}{0.9em}\\selectfont" << std::endl;
    out << "\\framebox{\\begin{tabular*}{60pt}{[email protected]{\\extracolsep{\\fill}}r}" << std::endl;
    
    // This is the particular info for the histogram
    out << "Entries & " ;
    formatNumber(histogram->GetEntries(), out) << " \\\\" << std::endl;
    out << "Mean & " ;
    formatNumber(histogram->GetMean(), out) << " \\\\" << std::endl;
    out << "RMS & " ;
    formatNumber(histogram->GetRMS(), out) << " \\\\" << std::endl;
    if (chi2 > 0.0) {
        out << "Fit $\\chi^{2}$ & " ;
        formatNumber(chi2, out) << " \\\\" << std::endl;
    }
    RooRealVar *theVariable;
    for (int index = 0; index < variables.getSize(); index++) {
        theVariable = dynamic_cast<RooRealVar*>(variables.find(variables[index].GetName()));
        out << theVariable->GetTitle() << " & $\\textrm{" ;
        formatNumber(theVariable->getValV(), out) << "} \\pm \\textrm{" ;
        formatNumber(theVariable->getError(), out) << "}$ \\\\" << std::endl;
    }
    out << "\\end{tabular*}}}" << std::endl;
    out << "\\end{document}" << std::endl;
    histogram->Delete();
}

void my_hook()
{
   hadaq::TrbProcessor* trb3 = base::ProcMgr::instance()->FindProc("TRB0");
   if (trb3==0) return;

   printf("Do extra work NUM %u\n", trb3->NumSubProc());

   for (unsigned ntdc=0xC000;ntdc<0xC005;ntdc++) {

      hadaq::TdcProcessor* tdc = trb3->GetTDC(ntdc);
      if (tdc==0) continue;

      TH1* hist = (TH1*) tdc->GetChannelRefHist(1);

      printf("  TDC%u mean:%5.2f rms:%5.2f\n", tdc->GetID(), hist->GetMean(), hist->GetRMS());

      tdc->ClearChannelRefHist(1);
   }

}

// Returns the RMS including 96% of the histogram entries, cutting the tails:
Double_t getRMS96(char* hs, double truncation=96.) {

  bool debug = false;

  TH1 *h = (TH1*)gDirectory->Get(hs);
  if( h == NULL ){ cout << hs << " does not exist\n"; return 0; }

  // Total entries:
  double integral = h->GetEntries();
  int maxbin = h->GetMaximumBin();
  if(debug) cout << "entries=" << integral << " maxbin=" << maxbin << endl;

  double subrange_integral = h->GetBinContent(maxbin);
  int bin = 0;
  while(subrange_integral < truncation/100*integral) {
    bin++;
    // Add one bin to the left:
    subrange_integral += h->GetBinContent(maxbin-bin);
    // Add one bin to the right:
    subrange_integral += h->GetBinContent(maxbin+bin);
    if(debug) cout << "subrange " << (maxbin-bin) << "-" << (maxbin+bin) << ": entries=" << subrange_integral << endl;
  }
  if(debug) cout << "subrange " << (maxbin-bin) << "-" << (maxbin+bin) << " now has " << subrange_integral << " entries, this is " << (100.0*subrange_integral)/integral << "%" << endl;

  // Correct by overshoot bin:
  subrange_integral -= h->GetBinContent(maxbin+bin);
  subrange_integral -= h->GetBinContent(maxbin-bin);
  bin--;

  int binlow = maxbin-bin;
  int binhigh = maxbin+bin;
  if(debug) cout << "subrange " << (maxbin-bin) << "-" << (maxbin+bin) << " now has " << subrange_integral << " entries, this is " << (100.0*subrange_integral)/integral << "%" << endl;

  h->GetXaxis()->SetRange(binlow,binhigh); //to restrict range to bins binlow to binhigh
  double rms96 = h->GetRMS(); //will return the RMS within the axis range

  return rms96;
}

void my_hook()
{
  hadaq::HldProcessor* hld = base::ProcMgr::instance()->FindProc("HLD");

  cout << Form("hook counter: %d", hook_counter) << endl;





  // Calibraton Time //  
  if (width_counter==0)
    {
      // int chId = 5;
      // gSystem->Exec(Form("echo -e 'asdf \tasdf sadf sadf as f' >> ~/git_cahitugur/trb3/measurements/ToT/pulseWidth_scan.txt"));

      // gSystem->Exec(Form("echo -e '%5.3f \t%d' >> ~/git_cahitugur/trb3/measurements/ToT/pulseWidth_scan.txt", width_counter*1.667, chId));


      if (hook_counter==15)
	{
	  width_counter++;
	  hook_counter = 0;
	  continue;
	}
      hook_counter++;
      continue;
    }
  ////////////////////////////////////

  // Pulser Settings //
  if (hook_counter==0) 
    {
      printf("Set pulser\n");
      gSystem->Exec(Form("perl ~/git_cahitugur/scripts/telnet_command_send.pl -pt 0x%x -pr 0x1fff -st 0x11", width_counter));
    }

  else if(hook_counter == 10)
    {
      width_counter++;
      hook_counter = 0;
      continue;
    }
  hook_counter++;
  //////////////////////////////////////

  // Histogram generation //
  int tdcmap[2] = { 0xC000, 0xC008 };

  if(hook_counter == 2)
    {
      cout << Form("Clear histograms\n") << endl;
      for (int cnt=0;cnt<2;cnt++) {
	
	hadaq::TdcProcessor* tdc = hld ? hld->FindTDC(tdcmap[cnt]) : 0;
	if (tdc==0) { printf("DID NOT FOUND TDC\n"); return; }
	for (int chId=1;chId<chNumber;chId=chId+1){
	  //      tdc_1->ClearChannelRefHist(chId);
	  tdc->ClearHist(chId, 5);
	}
      }
    }
  else if(hook_counter == 10)
    {
      gSystem->Exec(Form("echo -e '# pulseWidth(ns) \tTDC \tCHANNEL \tMEAN(ns) \tRMS(ps)' >> ~/git_cahitugur/trb3/measurements/ProgrammableOscillator/pulseWidth_scan.txt"));
      for (int cnt=0;cnt<2;cnt++) {
	
	hadaq::TdcProcessor* tdc = hld ? hld->FindTDC(tdcmap[cnt]) : 0;
	if (tdc==0) { printf("DID NOT FOUND TDC\n"); return; }
	
	for (int chId=1;chId<chNumber;chId=chId+1){
	  //  TH1* hist = (TH1*) tdc_1->GetChannelRefHist(chId); // argument is the channel number
	  TH1* tot = (TH1*) tdc->GetHist(chId,5); // arguments are the channel number and the ToT histogram number
	  gSystem->Exec(Form("echo -e '%5.3f \t\t%x \t%d \t\t%5.3f \t\t%5.4f' >> ~/git_cahitugur/trb3/measurements/ProgrammableOscillator/pulseWidth_scan.txt", width_counter*1.667, tdc->GetID(), chId, tot->GetMean(), tot->GetRMS()*1000 ));
	}
      }
    }
 }

//.........这里部分代码省略.........
            else if (htype == RarityType) 
               sig->SetTitle( Form("TMVA Rarity for classifier: %s", methodTitle.Data()) );
            else if (htype == CompareType) 
               sig->SetTitle( Form("TMVA overtraining check for classifier: %s", methodTitle.Data()) );
         
            // create new canvas
            TString ctitle = ((htype == MVAType) ? 
                              Form("TMVA response %s",methodTitle.Data()) : 
                              (htype == ProbaType) ? 
                              Form("TMVA probability %s",methodTitle.Data()) :
                              (htype == CompareType) ? 
                              Form("TMVA comparison %s",methodTitle.Data()) :
                              Form("TMVA Rarity %s",methodTitle.Data()));
         
            TString cname = ((htype == MVAType) ? 
                             Form("output_%s",methodTitle.Data()) : 
                             (htype == ProbaType) ? 
                             Form("probability_%s",methodTitle.Data()) :
                             (htype == CompareType) ? 
                             Form("comparison_%s",methodTitle.Data()) :
                             Form("rarity_%s",methodTitle.Data()));

            c = new TCanvas( Form("canvas%d", countCanvas+1), ctitle, 
                             countCanvas*50+200, countCanvas*20, width, (Int_t)width*0.78 ); 
    
            // set the histogram style
            TMVAGlob::SetSignalAndBackgroundStyle( sig, bgd );
   
            // normalise both signal and background
            TMVAGlob::NormalizeHists( sig, bgd );
   
            // frame limits (choose judicuous x range)
            Float_t nrms = 4;
            cout << "--- Mean and RMS (S): " << sig->GetMean() << ", " << sig->GetRMS() << endl;
            cout << "--- Mean and RMS (B): " << bgd->GetMean() << ", " << bgd->GetRMS() << endl;
            Float_t xmin = TMath::Max( TMath::Min(sig->GetMean() - nrms*sig->GetRMS(), 
                                                  bgd->GetMean() - nrms*bgd->GetRMS() ),
                                       sig->GetXaxis()->GetXmin() );
            Float_t xmax = TMath::Min( TMath::Max(sig->GetMean() + nrms*sig->GetRMS(), 
                                                  bgd->GetMean() + nrms*bgd->GetRMS() ),
                                       sig->GetXaxis()->GetXmax() );
            Float_t ymin = 0;
            Float_t maxMult = (htype == CompareType) ? 1.3 : 1.2;
            Float_t ymax = TMath::Max( sig->GetMaximum(), bgd->GetMaximum() )*maxMult;
   
            // build a frame
            Int_t nb = 500;
            TString hFrameName(TString("frame") + methodTitle);
            TObject *o = gROOT->FindObject(hFrameName);
            if(o) delete o;
            TH2F* frame = new TH2F( hFrameName, sig->GetTitle(), 
                                    nb, xmin, xmax, nb, ymin, ymax );
            frame->GetXaxis()->SetTitle( methodTitle + ((htype == MVAType || htype == CompareType) ? " response" : "") );
            if      (htype == ProbaType  ) frame->GetXaxis()->SetTitle( "Signal probability" );
            else if (htype == RarityType ) frame->GetXaxis()->SetTitle( "Signal rarity" );
            frame->GetYaxis()->SetTitle("Normalized");
            TMVAGlob::SetFrameStyle( frame );
   
            // eventually: draw the frame
            frame->Draw();  
    
            c->GetPad(0)->SetLeftMargin( 0.105 );
            frame->GetYaxis()->SetTitleOffset( 1.2 );

            // Draw legend               
            TLegend *legend= new TLegend( c->GetLeftMargin(), 1 - c->GetTopMargin() - 0.12, 

void makePlot(const std::string& inputFilePath, const std::string& canvasName, const std::string& sample, int massPoint, const std::string& channel, double k, 
	      const std::string& inputFileName, const std::string& outputFilePath, const std::string& outputFileName)
{
  std::string inputFileName_full = Form("%s%s", inputFilePath.data(), inputFileName.data());
  TFile* inputFile = new TFile(inputFileName_full.data());
  if ( !inputFile ) {
    std::cerr << "Failed to open input file = " << inputFileName_full << " !!" << std::endl;
    assert(0);
  }

  inputFile->ls();

  TCanvas* canvas = dynamic_cast<TCanvas*>(inputFile->Get(canvasName.data()));
  if ( !canvas ) {
    std::cerr << "Failed to load canvas = " << canvasName << " !!" << std::endl;
    assert(0);
  }

  int idxPad = -1;
  if ( massPoint ==  90 ) idxPad = 1;
  if ( massPoint == 125 ) idxPad = 2;
  if ( massPoint == 200 ) idxPad = 3;
  if ( massPoint == 300 ) idxPad = 4;
  if ( massPoint == 500 ) idxPad = 5;
  if ( massPoint == 800 ) idxPad = 6;  
  if ( !(idxPad >= 1 && idxPad <= 6) ) {
    std::cerr << "Invalid sample = " << sample << " !!" << std::endl;
    assert(0);
  }
  TVirtualPad* pad = canvas->GetPad(idxPad);
  std::cout << "pad = " << pad << ": ClassName = " << pad->ClassName() << std::endl;

  TCanvas* canvas_new = new TCanvas("canvas_new", "canvas_new", 900, 800);
  canvas_new->SetFillColor(10);
  canvas_new->SetBorderSize(2);
  canvas_new->SetTopMargin(0.065);
  canvas_new->SetLeftMargin(0.17);
  canvas_new->SetBottomMargin(0.165);
  canvas_new->SetRightMargin(0.015);
  canvas_new->SetLogx(true);
  canvas_new->SetLogy(true);
  canvas_new->Draw();
  canvas_new->cd();

  //TList* pad_primitives = canvas->GetListOfPrimitives();
  TList* pad_primitives = pad->GetListOfPrimitives();

  TH1* histogramCA            = 0;
  TH1* histogramSVfit         = 0;
  TH1* histogramSVfitMEMkEq0  = 0;
  TH1* histogramSVfitMEMkNeq0 = 0;

  TIter pad_nextObj(pad_primitives);
  while ( TObject* obj = pad_nextObj() ) {
    std::string objName = "";
    if ( dynamic_cast<TNamed*>(obj) ) objName = (dynamic_cast<TNamed*>(obj))->GetName();    
    std::cout << "obj = " << obj << ": name = " << objName << ", type = " << obj->ClassName() << std::endl;

    TH1* tmpHistogram = dynamic_cast<TH1*>(obj);
    if ( tmpHistogram ) {
      std::cout << "tmpHistogram:" 
		<< " fillColor = " << tmpHistogram->GetFillColor() << ", fillStyle = " << tmpHistogram->GetFillStyle() << ","
		<< " lineColor = " << tmpHistogram->GetLineColor() << ", lineStyle = " << tmpHistogram->GetLineStyle() << ", lineWidth = " << tmpHistogram->GetLineWidth() << ","
		<< " markerColor = " << tmpHistogram->GetMarkerColor() << ", markerStyle = " << tmpHistogram->GetMarkerStyle() << ", markerSize = " << tmpHistogram->GetMarkerSize() << ","
		<< " integral = " << tmpHistogram->Integral() << std::endl;
      std::cout << "(mean = " << tmpHistogram->GetMean() << ", rms = " << tmpHistogram->GetRMS() << ": rms/mean = " << (tmpHistogram->GetRMS()/tmpHistogram->GetMean()) << ")" << std::endl;
      if ( tmpHistogram->GetLineColor() == 416 ) histogramCA            = tmpHistogram;
      if ( tmpHistogram->GetLineColor() == 600 ) histogramSVfit         = tmpHistogram;
      if ( tmpHistogram->GetLineColor() == 616 ) histogramSVfitMEMkEq0  = tmpHistogram;
      if ( tmpHistogram->GetLineColor() == 632 ) histogramSVfitMEMkNeq0 = tmpHistogram;
    }
  }

  if ( !(histogramCA && histogramSVfit && histogramSVfitMEMkEq0 && histogramSVfitMEMkNeq0) ) {
    std::cerr << "Failed to load histograms !!" << std::endl;
    assert(0);
  }

  //gStyle->SetLineStyleString(2,"40 10 10 10 10 10 10 10");
  //gStyle->SetLineStyleString(3,"25 15");
  //gStyle->SetLineStyleString(4,"60 25");

  //int colors[4] = { kBlack, kGreen - 6, kBlue - 7, kMagenta - 7  };
  int colors[4] = { 28, kGreen - 6, kBlue - 7, kBlack };
  //int lineStyles[4] = { 2, 3, 4, 1 };
  int lineStyles[4] = { 7, 1, 1, 1 };
  //int lineWidths[4] = { 3, 3, 4, 3 };
  int lineWidths[4] = { 3, 3, 1, 1 };
  int markerStyles[4] = { 20, 25, 21, 24 };
  int markerSizes[4] = { 2, 2, 2, 2 };

  histogramCA->SetFillColor(0);
  histogramCA->SetFillStyle(0);
  histogramCA->SetLineColor(colors[0]);
  histogramCA->SetLineStyle(lineStyles[0]);
  histogramCA->SetLineWidth(lineWidths[0]);
  histogramCA->SetMarkerColor(colors[0]);
  histogramCA->SetMarkerStyle(markerStyles[0]);
  histogramCA->SetMarkerSize(markerSizes[0]);

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

void printCalibStat(Int_t run, const char * fname,  TTreeSRedirector * pcstream){

  //
  // Dump the statistical information about all histograms in the calibration files 
  //    into the statistical tree, print on the screen (log files) as well 
  //
  //
  // 1. Default dump for all histograms
  //    Information to dump:
  //    stat =Entries, Mean, MeanError,  RMS, MaxBin
  //    Branch naming convention:
  //    <detName>_<hisName><statName>
  //
  // 2. Detector statistical information  - to be implemented by expert
  //                                      - First version implemented by MI 
  //  
  // 

  TFile *fin = TFile::Open(fname);
  if (!fin) return;
  const Double_t kMaxHis=10000;
  
  TList * keyList = fin->GetListOfKeys();
  Int_t nkeys=keyList->GetEntries();
  Double_t *hisEntries = new Double_t[kMaxHis];
  Double_t *hisMean = new Double_t[kMaxHis];
  Double_t *hisMeanError = new Double_t[kMaxHis];
  Double_t *hisRMS = new Double_t[kMaxHis];
  Double_t *hisMaxBin = new Double_t[kMaxHis];
  Int_t counter=0;
  
  if (pcstream) (*pcstream)<<"calibStatAll"<<"run="<<run;
  for (Int_t ikey=0; ikey<nkeys; ikey++){
    TObject * object = fin->Get(keyList->At(ikey)->GetName());
    if (!object) continue;
    if (object->InheritsFrom("TCollection")==0) continue;
    TSeqCollection *collection  = (TSeqCollection*)object; 
    Int_t nentries= collection->GetEntries();
    for (Int_t ihis=0; ihis<nentries; ihis++){
      TObject * ohis = collection->At(ihis);
      if (!ohis) continue;
      if (ohis->InheritsFrom("TH1")==0) continue;
      TH1* phis = (TH1*)ohis;
      hisEntries[counter]=phis->GetEntries();	
      Int_t idim=1;
      if (ohis->InheritsFrom("TH2")) idim=2;
      if (ohis->InheritsFrom("TH3")) idim=3;        
      hisMean[counter]=phis->GetMean(idim);	
      hisMeanError[counter]=phis->GetMeanError(idim);	
      hisRMS[counter]=phis->GetRMS(idim);	
      hisMaxBin[counter]=phis->GetBinCenter(phis->GetMaximumBin());	
      if (pcstream) (*pcstream)<<"calibStatAll"<<
		      Form("%s_%sEntries=",keyList->At(ikey)->GetName(), phis->GetName())<<hisEntries[counter]<<	
		      Form("%s_%sMean=",keyList->At(ikey)->GetName(), phis->GetName())<<hisMean[counter]<<	
		      Form("%s_%sMeanError=",keyList->At(ikey)->GetName(), phis->GetName())<<hisMeanError[counter]<<	
		      Form("%s_%sRMS=",keyList->At(ikey)->GetName(), phis->GetName())<<hisRMS[counter]<<	
		      Form("%s_%sMaxBin=",keyList->At(ikey)->GetName(), phis->GetName())<<hisMaxBin[counter];	
      //printf("Histo:\t%s_%s\t%f\t%d\n",keyList->At(ikey)->GetName(), phis->GetName(), hisEntries[counter],idim);
      counter++;
    }
    delete object;
  }    
  
  //
  // Expert dump example (MI first iteration):
  //
  // 0.)  TOF dump
  //

  Int_t tofEvents=0;
  Int_t tofTracks=0;
  TList * TOFCalib = (TList*)fin->Get("TOFHistos");      
  if (TOFCalib) {
    TH1 *histoEvents = (TH1*)TOFCalib->FindObject("hHistoVertexTimestamp");
    TH1 *histoTracks = (TH1*)TOFCalib->FindObject("hHistoDeltatTimestamp");
    if (histoEvents && histoTracks){
      tofEvents = TMath::Nint(histoEvents->GetEntries());
      tofTracks = TMath::Nint(histoTracks->GetEntries());
    }
    delete TOFCalib;
  }
  printf("Monalisa TOFevents\t%d\n",tofEvents);
  if (pcstream) (*pcstream)<<"calibStatAll"<<"TOFevents="<<tofEvents;
  printf("Monalisa TOFtracks\t%d\n",tofTracks);
  if (pcstream) (*pcstream)<<"calibStatAll"<<"TOFtracks="<<tofTracks;

  //
  // 1.)  TPC  dump - usefull events/tracks  for the calibration
  //
  Int_t tpcEvents=0;
  Int_t tpcTracks=0;
  TObject* obj = dynamic_cast<TObject*>(fin->Get("TPCCalib"));
  TObjArray* array = dynamic_cast<TObjArray*>(obj);
  TDirectory* dir = dynamic_cast<TDirectory*>(obj);
  AliTPCcalibTime  * calibTime = NULL;
  if (dir) {
    calibTime = dynamic_cast<AliTPCcalibTime*>(dir->Get("calibTime"));
  }
  else if (array){
    calibTime = (AliTPCcalibTime *)array->FindObject("calibTime");
//.........这里部分代码省略.........

TGraph* makeGraph(std::vector<histogram_vs_X_Type>& histograms_vs_X, const std::vector<double>& y_true_array, const std::string& mode)
{
  //std::cout << "<makeGraph>:" << std::endl;

  enum { kResponse, kResponse_asymmetric, kResolution, kResolution_asymmetric_left, kResolution_asymmetric_right, kFitResponse, kFitResolution };
  int mode_int = -1;
  if      ( mode == "response"                    ) mode_int = kResponse;
  else if ( mode == "response_asymmetric"         ) mode_int = kResponse_asymmetric;
  else if ( mode == "resolution"                  ) mode_int = kResolution;
  else if ( mode == "resolution_asymmetric_left"  ) mode_int = kResolution_asymmetric_left;
  else if ( mode == "resolution_asymmetric_right" ) mode_int = kResolution_asymmetric_right;
  else if ( mode == "fit_resolution"              ) mode_int = kFitResolution;
  else if ( mode == "fit_resolution"              ) mode_int = kFitResolution;
  else assert(0);

  unsigned numPoints = histograms_vs_X.size();
  //std::cout << " numPoints = " << numPoints << std::endl;
  assert(numPoints > 0);

  TGraphAsymmErrors* graph = new TGraphAsymmErrors(numPoints);

  for ( unsigned iPoint = 0; iPoint < numPoints; ++iPoint ) {
    double x = histograms_vs_X[iPoint].x_;
    double xErrUp = histograms_vs_X[iPoint].xErrUp_;
    double xErrDown = histograms_vs_X[iPoint].xErrDown_;

    TH1* histogram = histograms_vs_X[iPoint].histogram_;

    double histogram_mean = histogram->GetMean();
    double histogram_meanErr = histogram->GetMeanError();
    double histogram_rms = histogram->GetRMS();
    double histogram_rmsErr = histogram->GetRMSError();

    double y_true = y_true_array[iPoint];

    double y, yErrUp, yErrDown;
    if ( mode_int == kResponse ) {
      y = histogram_mean/y_true;
      yErrUp = histogram_meanErr/y_true;
      yErrDown = yErrUp;
    } else if ( mode_int == kResponse_asymmetric ) {
      return makeGraph_response_asymmetric(histograms_vs_X, y_true_array);
    } else if ( mode_int == kResolution ) {
      y = histogram_rms/histogram_mean;
      yErrUp = y*TMath::Sqrt(TMath::Power(histogram_rms/histogram_rmsErr, 2.) + TMath::Power(histogram_mean/histogram_meanErr, 2.));
      yErrDown = yErrUp;
    } else if ( mode_int == kResolution_asymmetric_left ) {
      return makeGraph_resolution_asymmetric(histograms_vs_X).first;
    } else if ( mode_int == kResolution_asymmetric_right ) {
      return makeGraph_resolution_asymmetric(histograms_vs_X).second;
    } else {
      std::vector<fitParameterType> fitParameter;
      fit(histogram, -1., -1., y_true, fitParameter);
      fitParameterType fitParameterOfInterest;
      if      ( mode_int == kFitResolution ) fitParameterOfInterest = fitParameter[2];
      else if ( mode_int == kFitResponse   ) fitParameterOfInterest = fitParameter[3];
      else assert(0);
      y = fitParameterOfInterest.value_;
      yErrUp = fitParameterOfInterest.errUp_;
      yErrDown = fitParameterOfInterest.errDown_;
    }

    graph->SetPoint(iPoint, x, y);
    graph->SetPointError(iPoint, xErrDown, xErrUp, yErrDown, yErrUp);
  }

  return graph;
}

void EMCDistribution_SUM_RawADC(TString sTOWER = "Energy_Sum_col1_row2_5x5",
                                TString CherenkovSignal = "C2_Inner")
{
  TH1 *EnergySum_LG_full = new TH1F("EnergySum_LG_full",
                                    ";Tower Energy Sum (ADC);Count / bin", 260, -100, 2500);
  TH1 *EnergySum_LG = new TH1F("EnergySum_LG",
                               ";Tower Energy Sum (ADC);Count / bin", 260, -100, 2500);
  //  TH1 * EnergySum_HG = new TH1F("EnergySum_HG",
  //      ";Low range Tower Energy Sum (ADC);Count / bin", 50, 0, 500);

  TH1 *C2_Inner_full = new TH1F("C2_Inner_full",
                                CherenkovSignal + ";Cherenkov Signal (ADC);Count / bin", 1000, 0, 2000);
  TH1 *C2_Inner = new TH1F("C2_Inner",
                           CherenkovSignal + ";Cherenkov Inner Signal (ADC);Count / bin", 1000, 0, 2000);

  EnergySum_LG_full->SetLineColor(kBlue + 3);
  EnergySum_LG_full->SetLineWidth(2);

  EnergySum_LG->SetLineColor(kGreen + 3);
  EnergySum_LG->SetLineWidth(3);
  EnergySum_LG->SetMarkerColor(kGreen + 3);

  C2_Inner_full->SetLineColor(kBlue + 3);
  C2_Inner_full->SetLineWidth(2);

  C2_Inner->SetLineColor(kGreen + 3);
  C2_Inner->SetLineWidth(3);
  C2_Inner->SetMarkerColor(kGreen + 3);

  TCut c2 = CherenkovSignal + ">240";

  T->Draw(sTOWER + ">>EnergySum_LG_full", "", "goff");
  T->Draw(sTOWER + ">>EnergySum_LG", c2, "goff");
  T->Draw(CherenkovSignal + ">>C2_Inner_full", "", "goff");
  T->Draw(CherenkovSignal + ">>C2_Inner", c2, "goff");

  TText *t;
  TCanvas *c1 = new TCanvas(
      "EMCDistribution_SUM_RawADC_" + sTOWER + "_" + CherenkovSignal + cuts,
      "EMCDistribution_SUM_RawADC_" + sTOWER + "_" + CherenkovSignal + cuts, 1800,
      600);
  c1->Divide(3, 1);
  int idx = 1;
  TPad *p;

  p = (TPad *) c1->cd(idx++);
  c1->Update();
  p->SetLogy();
  p->SetGridx(0);
  p->SetGridy(0);

  C2_Inner_full->DrawClone();
  C2_Inner->DrawClone("same");

  p = (TPad *) c1->cd(idx++);
  c1->Update();
  p->SetLogy();
  p->SetGridx(0);
  p->SetGridy(0);

  TH1 *h = (TH1 *) EnergySum_LG_full->DrawClone();
  //  h->GetXaxis()->SetRangeUser(0, h->GetMean() + 5 * h->GetRMS());
  (TH1 *) EnergySum_LG->DrawClone("same");

  p = (TPad *) c1->cd(idx++);
  c1->Update();
  //  p->SetLogy();
  p->SetGridx(0);
  p->SetGridy(0);

  TH1 *h_full = (TH1 *) EnergySum_LG_full->DrawClone();
  TH1 *h = (TH1 *) EnergySum_LG->DrawClone("same");

  TF1 *fgaus_g = new TF1("fgaus_LG_g", "gaus", h->GetMean() - 1 * h->GetRMS(),
                         h->GetMean() + 4 * h->GetRMS());
  fgaus_g->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
                         h->GetMean() + 2 * h->GetRMS());
  h->Fit(fgaus_g, "MR0N");

  TF1 *fgaus = new TF1("fgaus_LG", "gaus",
                       fgaus_g->GetParameter(1) - 1 * fgaus_g->GetParameter(2),
                       fgaus_g->GetParameter(1) + 4 * fgaus_g->GetParameter(2));
  fgaus->SetParameters(fgaus_g->GetParameter(0), fgaus_g->GetParameter(1),
                       fgaus_g->GetParameter(2));
  h->Fit(fgaus, "MR");

  h->Sumw2();
  h_full->Sumw2();
  h_full->GetXaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
                                   h->GetMean() + 4 * h->GetRMS());

  h->SetLineWidth(2);
  h->SetMarkerStyle(kFullCircle);

  h_full->SetTitle(
      Form("#DeltaE/<E> = %.1f%%",
           100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));

  //  p = (TPad *) c1->cd(idx++);
  //  c1->Update();
//.........这里部分代码省略.........

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

		for (int j = 0; j < 1; j++) {
			fread((void *) &event_data, 1, sizeof(event_data), fdata);
			p = (struct channel_struct *) &event_data.ch[0]; // read bunch of data
			p += anaChannel;
			for (ch = 0; ch < DOMINO_NCELL; ch++) {
				grPedData->GetPoint(ch, itmp, PedVal);
				reftmp = TMath::Abs((Double_t)(p->data[ch])-PedVal);
				if (reftmp > 0.8* refval)
					refval = 0.2*reftmp+0.8*refval;
				//					cout << ch << " " <<p->data[ch] << " " <<reftmp << " " << refval <<endl ;
			}
		}
		cout << "refval="<< refval<<endl;
		rewind(fdata);

		Int_t ievt = 1;
		// go to first event (startEv)
		while (ievt < startEv) {
			fread((void *) &event_data, 1, sizeof(event_data), fdata);
			if (feof(fdata))
				break;
			ievt++;
		}

		ievt = 1;
		Int_t iTrig = 0;
		Int_t flagEnd = 0;
		Double_t chtmp, chtrig;
		Double_t ratio;
		Double_t mean, rms;

		// loop on events

		while (ievt <= nevt && !flagEnd) {

			fread((void *) &event_data, 1, sizeof(event_data), fdata);
			if (feof(fdata))
				flagEnd = 1;

			p = (struct channel_struct *) &event_data.ch[0]; // read bunch of data
			dep = (struct channel_struct *) &event_data.ch[1]; // read bunch of data

			//now anaChannel analysis

			p += anaChannel;

			// read data, subtract pedestals values and fill hCellCalibList.

			for (int ch = 0; ch < DOMINO_NCELL; ch++) {
				// Read pedestal value for this cell
				grPedData->GetPoint(ch, itmp, PedVal);
				chtmp = (Double_t)(p->data[ch]); // data value
				chtmp = chtmp - PedVal;
				//if (TMath::Abs(chtmp) > 0.9 * refval)
				((TH1 *) hCellCalibList->At(iFile*DOMINO_NCELL+ch))->Fill(chtmp);
				//cout << ch << " " << iFile << " " << chtmp << endl;
			}

			gProgress->Increment(1);
			gSystem->DispatchOneEvent(kTRUE);
			ievt++; // next event
		}


		TH1 *hCellTmp;
		for(ch = 0; ch < DOMINO_NCELL;ch++) {
			hCellTmp = ((TH1 *) hCellCalibList->At(iFile*DOMINO_NCELL+ch));
			//hCellTmp->Fit("gaus", "Q");
			//mean = (hCellTmp->GetFunction("gaus"))->GetParameter(1);
			//rms = (hCellTmp->GetFunction("gaus"))->GetParameter(2);
			mean = hCellTmp->GetMean();
			rms = hCellTmp->GetRMS();
			((TGraphErrors *) (grCellCalibList->At(ch)))->SetPoint(iFile, (Double_t) mV[iFile], mean);
			((TGraphErrors *) (grCellCalibList->At(ch)))->SetPointError(iFile, 0., rms);
		}


		ctest->cd(iFile + 1);
		hCellTmp = ((TH1 *) hCellCalibList->At(567 + iFile*DOMINO_NCELL));
		hCellTmp->Fit("gaus","Q");
		hCellTmp->DrawCopy();
	}

	TString OutFile = "CalibrationDataNew";
	OutFile += nevt;
	OutFile += "events.root";
	TFile *f = new TFile(OutFile, "RECREATE");
	for (int ch = 0; ch < DOMINO_NCELL; ch++) {
		TString key = "CalibDataCell";
		key += ch;
		((TGraphErrors*) grCellCalibList->At(ch))->Write(key);
	}
	f->Close();

	hCellCalibList->Delete();

	((TGMainFrame *) gProgress->GetParent())->CloseWindow();
	fclose(fdata);
}

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

	gProgress->Reset();
	gProgress->SetMax(nevt);

	gSystem->ProcessEvents();

	while (ievt <= nevt && !flagEnd) {
		fread((void *) &event_data, 1, sizeof(event_data), fdata);
		if (feof(fdata))
			flagEnd = 1;

		p = (struct channel_struct *) &event_data.ch[0]; // read bunch of data
		dep = (struct channel_struct *) &event_data.ch[1]; // read bunch of data
		// goes to channel to analyze
		p += anaChannel;

		// read data, subtract pedestals values and save results in grAnaChDataTemp graph with
		// fixed error for each point (x = 0.5 and y = 2.1). Also generate an array with Domino
		// X and Y values

		TGraphErrors *grAnaChDataTemp = new TGraphErrors(DOMINO_NCELL);

		for (int ch = 0; ch < DOMINO_NCELL; ch++) {
			// Read pedestal value for this cell
			grPed->GetPoint(ch, itmp, PedVal);
			chtmp = (Double_t)(p->data[ch]); // data value
			chtmp = chtmp - PedVal;
			grAnaChDataTemp->SetPoint(ch, (Double_t) ch, chtmp);
			grAnaChDataTemp->SetPointError(ch, 0.5, 2.1);
		}
		// create fit functions
		TF1 *fsin = new TF1("fsin", sigSin, 0., 1024., 4);
		fsin->SetParameters(600., 255., 150., 150.);
		fsin->SetParNames("amplitude", "Period", "Phase", "DC-Offset");

		grAnaChDataTemp->Fit("fsin", "Q");
		TF1 *fsinFit = grAnaChDataTemp->GetFunction("fsin");
		fsinFit->SetParNames("amplitude", "Period", "Phase", "DC-Offset");

		// debug
		cfitTest->cd(ievt);
		grAnaChDataTemp->SetMarkerStyle(20);
		grAnaChDataTemp->SetMarkerSize(0.3);
		grAnaChDataTemp->GetYaxis()->SetLabelSize(0.12);
		grAnaChDataTemp->GetXaxis()->SetLabelSize(0.12);
		grAnaChDataTemp->Draw("APE");

		Double_t fitPeriod, fitAmplitude, chisquare;
		fitPeriod = fsinFit->GetParameter("Period");
		fitAmplitude = TMath::Abs(fsinFit->GetParameter("amplitude"));
		chisquare = fsinFit->GetChisquare();

		cout << "period: " << fitPeriod << " amplitude: " << fitAmplitude << " chisquare: " << chisquare << endl;

		if(chisquare > 0.1e+06) {
			gProgress->Increment(1);
			gSystem->DispatchOneEvent(kTRUE);
			ievt++;
			continue;
		}

		gProgress->Increment(1);
		gSystem->DispatchOneEvent(kTRUE);

		hPeriod->Fill(fitPeriod);
		fitusati++;

		ievt++;

	}

	cout << "fit scartati :" << nevt - fitusati << endl;
	//draw
	TString Title = "Period distribution for nevt events";
	TCanvas *cPeriod = new TCanvas("cPeriod", Title, 700, 700);
	hPeriod->Draw();
	hPeriod->Fit("gaus");

	TF1 *fgausFit = hPeriod->GetFunction("gaus");
	//mean = fgausFit->GetParameter(1);
	//	rms = fgausFit->GetParameter(2);

	mean = hPeriod->GetMean();
	rms = hPeriod->GetRMS();

	TString OutFile = "Period";
	OutFile += nevt;
	OutFile += "events.dat";
	FILE *f = fopen(OutFile.Data(), "w");
	fwrite(&mean, sizeof(mean), 1, f);
	fwrite(&rms, sizeof(rms), 1, f);

	((TGMainFrame *) gProgress->GetParent())->CloseWindow();
	fclose(f);

	cout << "mean: " << mean << " rms: " << rms << endl;

	fclose(fdata);

}