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

//------------------------------------------------------------------------
void DrawQTC()
{
  
  TCanvas *c1 = new TCanvas("c1", "QTC C side",0,48,1280,951);
   c1->Divide(4,3);
   // c1->Divide(2,3);
  
  Char_t buf1[10];
  for (Int_t i=0; i<12; i++)
    {
      c1->cd(i+1);
      sprintf(buf1,"QTC%i",i+1);
      
      TH1F *qtc = (TH1F*) gFile->Get(buf1);
      
      Float_t mean = qtc->GetMean();
      Float_t rms = qtc->GetRMS();
      Float_t hminR=mean - 0.1*mean;
      Float_t hmaxR =mean + 0.1*mean;
      qtc->GetXaxis()->SetRange(hminR,hmaxR);
      Float_t hmin=mean - 3*rms;
      Float_t hmax =mean + 3*rms;
      qtc->GetXaxis()->SetRange(hmin,hmax);

     // TF1 *g2 = new TF1("g2", "gaus", hmin, hmax);
       //             qtc->Fit("g2","RQ");
	qtc->GetXaxis()->SetLabelSize(0.03);
       
       qtc->Draw();
    }
  
  TCanvas *c2 = new TCanvas("c2", "QTC A side",0,48,1280,951);
   c2->Divide(4,3);
   // c1->Divide(2,3);
  
  Char_t buf1[10];
  for (Int_t i=12; i<24; i++)
    {
      c2->cd(i+1-12);
      sprintf(buf1,"QTC%i",i+1);
      
      TH1F *qtc = (TH1F*) gFile->Get(buf1);
      
      Float_t mean = qtc->GetMean();
      Float_t rms = qtc->GetRMS();
      Float_t hminR=mean - 0.1*mean;
      Float_t hmaxR =mean + 0.1*mean;
      qtc->GetXaxis()->SetRange(hminR,hmaxR);
      Float_t hmin=mean - 3*rms;
      Float_t hmax =mean + 3*rms;
      qtc->GetXaxis()->SetRange(hmin,hmax);
      // TF1 *g2 = new TF1("g2", "gaus", hmin, hmax);
       //             qtc->Fit("g2","RQ");
 	qtc->GetXaxis()->SetLabelSize(0.03);
      
       qtc->Draw();
    }
  
  
}

vector<double> one_fit(bool do_fit)
{
    TH1F* h = new TH1F("h","h", 20, -10, 10);
    h->FillRandom("gaus");
    vector<double> ret;
    if (do_fit) {
        h->Fit("gaus","");
        TF1* g = h->GetFunction("gaus");
        h->Draw();
        canvas->Modified();
        canvas->Update();
        for (int ind=0; ind < 3; ++ind) {
            ret.push_back(g->GetParameter(ind));
        }
    }
    else {
        for (int ind=0; ind < 3; ++ind) {
            ret.push_back(0.0);
        }
    }
    
    ret.push_back(h->GetMean());
    ret.push_back(h->GetRMS());
    delete h;
    return ret;
}

double Chip::getMeanValueOfStrips(){
  double value=0;
  double maxValue=0;
  //TH1F * Histo = new TH1F("name","title",1000,0,);
  
  for(int i = 1 ; i <= 8 ; i++){
    if(this->getChipStrip(i)->getRate() > maxValue){
      maxValue = this->getChipStrip(i)->getRate();
    }
  }
  
  TH1F * Histo = new TH1F("name","title",1000,0,maxValue);
  for(int i = 1 ; i <= 8;i++){
    if( ! this->getChipStrip(i)->getIsDisconnected()){
      Histo->Fill(this->getChipStrip(i)->getRate());
    }
  }
  
  value = Histo->GetMean();
  delete Histo;
  if (maxValue == 0){
    cout << "maxValue is zero ! " << endl;
    value = maxValue;
  }
  
  return value;
}

void KVElasticCountRates::PrintResults(Double_t beam_intensity)
{
   // Print mean energy deposit & counting rate for given beam intensity in particles per second

   TIter it(&fHistos);
   TH1F* h;
   fRates.clear();

   std::vector<count_rate> count_rates;

   while ((h = (TH1F*)it())) {
      TString name = h->GetName();
      if (!name.EndsWith("_dW") && !name.EndsWith("_map")) {
         TH2F* map = (TH2F*)fHistos.FindObject(name + "_map");
         double rate = h->Integral() * fAtomicDensity * beam_intensity * fVolume / fNtirages;
         double emean = h->GetMean();
         KVDetector* det = gMultiDetArray->GetDetector(name);
         double fluence = rate / det->GetEntranceWindowSurfaceArea();
         double dissipation = emean * rate / det->GetEntranceWindowSurfaceArea();
         count_rates.push_back(
            count_rate(name, rate, emean, map->GetMean(), map->GetMean(2), fluence, dissipation)
         );
         fRates[name.Data()] = KVElasticCountRate(rate, emean, fluence, dissipation);
      }
   }
   std::sort(count_rates.begin(), count_rates.end(), compare_count_rates);

   for (std::vector<count_rate>::iterator it = count_rates.begin(); it != count_rates.end(); ++it) {
      it->print();
   }
}

void explore_zmass_boost()
{
  // Tell root not to draw everything to the screen.
  gROOT->SetBatch();

  //    TString dyfilename         = TString("/home/acarnes/h2mumu/samples/stage1/monte_carlo/bg/stage_1_dy_jetsToLL_asympt50_ALL.root");
  TString dyfilename         = TString("/home/acarnes/h2mumu/samples/stage1/monte_carlo/bg/stage_1_dy_ZToMuMu_asympt50_ALL.root");
  TString datafilename         = TString("/home/acarnes/h2mumu/samples/stage1/data_from_json/Cert_246908-251883_13TeV_PromptReco_Collisions15_JSON_v2/stage_1_doubleMuon_RunBPrompt_MINIAOD.root");

  TString savedir = TString("../png/dy_vs_data/run1cuts_v2_golden_json/dyzmumu/");

  // Initialize the DiMuPlottingSystems for MC and data
  DiMuPlottingSystem* dpsdata = new DiMuPlottingSystem(datafilename);
  DiMuPlottingSystem* dpsdy = new DiMuPlottingSystem(dyfilename);
  addDiMuMassPrimeBranch(dpsdata);
  dpsdata->applyRun1Cuts();
  dpsdy->applyRun1Cuts();

  // Get the 2D histos
  TH2F* hdataZPt = ZMassVsZPtHist2D("data_zpt", "recoCandMassPrime", "(14,0,60,30,87,95)", dpsdata);
  TH2F* hdyZPt = ZMassVsZPtHist2D("dy_zpt", "recoCandMass", "(14,0,60,30,87,95)", dpsdy);

  //overlayTProfiles(TH2F* hdata, TH2F* hdy, TString bininfo, TString savename)
  TProfile* p = overlayTProfiles(hdataZPt, hdyZPt, 90, 92, savedir+"z_mass_vs_z_pt.png");
  gStyle->SetOptFit(0011);
  fitTProfileCustom(p);
  TCanvas* c = new TCanvas();
  c->cd();
  p->Draw("hist c");
  p->Draw("E0 X0 same");
  dpsdata->arrangeStatBox(c);
  c->Draw();
  c->Print("blah2.png");

  // Look at RMS
  c->Clear();
  c->cd();
  c->SetGridx(kTRUE);
  c->SetGridy(kTRUE);
  c->cd();
  TH1F* h = dpsdata->hist1D("recoCandMass", "(100,87,95)", "");
  TH1F* h2 = dpsdata->hist1D("recoCandMassPrime", "(100,87,95)", "");
  h2->SetLineColor(2);
  h->Draw("");
  h2->Draw("same");

  std::cout << "mass  rms: " << h->GetRMS() << std::endl;
  std::cout << "mass' rms: " << h2->GetRMS() << std::endl;
  std::cout << "mass  mean: " << h->GetMean() << std::endl;
  std::cout << "mass' mean: " << h2->GetMean() << std::endl;
  std::cout << "mass  num: " << h->Integral() << std::endl;
  std::cout << "mass' num: " << h2->Integral() << std::endl;

  DiMuPlottingSystem* dps = new DiMuPlottingSystem();
  dps->arrangeStatBox(c);
  c->Draw();
  c->Print("newfit.png");
}

void PlotHistsNhitsPerModule(TFile* f, TTree* tr, TString strMillepedeRes, TString strOutdir)
{
  TString canvName="c_";
  canvName+=strMillepedeRes;
  canvName+="_";
  canvName+=StrPlotType(NHITS);
  canvName.ReplaceAll(".res","");


  //enum {PXB,PXF,TIB,TID,TOB,TEC};
  int colors[6]={1,2,3,4,6,7};
//  TString labels[6]={"PXB","PXF","TIB","TID","TOB","TEC"};

  f->cd();
  TCanvas* canv = new TCanvas(canvName,canvName,600,600);
  canv->SetLogx();
  canv->SetLogy();

  for (int ind=1; ind<=1; ind++){
    TString strHist = "hNhits_";
    strHist+=StrPar(ind);
    TString strCut="label<700000 && ((label%20-1)%9+1)==";
    strCut+=ind;
    TStyle style; 
    style.SetTitleFontSize(0.2);
    THStack *hSt = new THStack("hNhits","# of derivatives (~tracks or hits) per module");
    TLegend *leg = new TLegend(0.75,0.65,0.95,0.95);
    for (int inv=0; inv<6; inv++){
      std::cout<<"- - - - - -"<<std::endl;
      std::cout<<subdLabels[inv]<<":"<<std::endl;
      std::cout<<StrCutSubd(inv)<<": "<<tr->GetEntries(StrCutSubd(inv))<<" parameters"<<std::endl;
      TString strHist1=strHist;
      strHist1+=ind;
      strHist1+=inv;
      TH1F* hValInt = new TH1F(strHist1,strHist1,300,10,15000);  
      TString strCut1 = strCut+TString(" && ")+StrCutSubd(inv);
      tr->Draw(TString("Nhits>>")+strHist1,strCut1,"goff");
      std::cout<<"# hits = "<<(int)hValInt->GetMean()<<"+-"<<(int)hValInt->GetRMS()<<std::endl;
      hValInt->SetLineColor(1);
      hValInt->SetFillColor(colors[inv]);
      hValInt->SetLineWidth(2);
      hSt->Add(hValInt);
      leg->AddEntry(hValInt,subdLabels[inv],"f");
      leg->SetFillColor(0);
    }
    hSt->Draw();
    leg->Draw("same");
    
  }//end of loop over ind

  canvName+=".png";
  TString saveName=strOutdir+canvName;
  canv->SaveAs(saveName);
  saveName.ReplaceAll(".png",".pdf");
  canv->SaveAs(saveName);
}//end of PlotHistsNhitsPerModule

double GetVariable(const char *var,TCut cut)
{
  TTree *t=(TTree*) gROOT->FindObject("t");

  TH1F *htemp = (TH1F*) gROOT->FindObject("htemp");
  if(htemp) delete htemp;
  t->Draw(Form("%s>>htemp",var),cut);
  htemp = (TH1F*) gROOT->FindObject("htemp");
  return double(htemp->GetMean());
}

void GetChamberIDs(int IDArray[9]){
TCanvas *IDcanv = new TCanvas ("idGraph", "idGraph");
IDcanv->cd();
TH1F *idDummy = new TH1F("idDummy", "idDummy", 10, 220000000, 221000000);
idDummy->Draw();
for (int chamber=0; chamber<9; ++chamber){
  TString idCut = Form ("cham==%d", chamber);
  Calibration->Project("idDummy", "id", idCut);
  Int_t idNum = idDummy->GetMean();
  IDArray[chamber]=idNum;
}
}

void billtr(Int_t compress) {
   //read N histograms from a tree
   timer.Start();
   TFile f("billt.root");
   TH1F *h = 0;
   TTree *T = (TTree*)f.Get("T");
   T->SetBranchAddress("event",&h);
   TH1F *hmeant = new TH1F("hmeant","hist mean from tree",100,0,1);
   Long64_t nentries = T->GetEntries();
   for (Long64_t i=0;i<nentries;i++) {
      T->GetEntry(i);
      hmeant->Fill(h->GetMean());
   }
   timer.Stop();
   printf("billtr%d : RT=%7.3f s, Cpu=%7.3f s\n",compress,timer.RealTime(),timer.CpuTime());
}

void langaus() {
   // Fill Histogram
   Int_t data[100] = {0,0,0,0,0,0,2,6,11,18,18,55,90,141,255,323,454,563,681,
                    737,821,796,832,720,637,558,519,460,357,291,279,241,212,
                    153,164,139,106,95,91,76,80,80,59,58,51,30,49,23,35,28,23,
                    22,27,27,24,20,16,17,14,20,12,12,13,10,17,7,6,12,6,12,4,
                    9,9,10,3,4,5,2,4,1,5,5,1,7,1,6,3,3,3,4,5,4,4,2,2,7,2,4};
   TH1F *hSNR = new TH1F("snr","Signal-to-noise",400,0,400);

   for (Int_t i=0; i<100; i++) hSNR->Fill(i,data[i]);

   // Fitting SNR histo
   printf("Fitting...\n");

   // Setting fit range and start values
   Double_t fr[2];
   Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
   fr[0]=0.3*hSNR->GetMean();
   fr[1]=3.0*hSNR->GetMean();

   pllo[0]=0.5; pllo[1]=5.0; pllo[2]=1.0; pllo[3]=0.4;
   plhi[0]=5.0; plhi[1]=50.0; plhi[2]=1000000.0; plhi[3]=5.0;
   sv[0]=1.8; sv[1]=20.0; sv[2]=50000.0; sv[3]=3.0;

   Double_t chisqr;
   Int_t    ndf;
   TF1 *fitsnr = langaufit(hSNR,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);

   Double_t SNRPeak, SNRFWHM;
   langaupro(fp,SNRPeak,SNRFWHM);

   printf("Fitting done\nPlotting results...\n");

   // Global style settings
   gStyle->SetOptStat(1111);
   gStyle->SetOptFit(111);
   gStyle->SetLabelSize(0.03,"x");
   gStyle->SetLabelSize(0.03,"y");

   hSNR->GetXaxis()->SetRange(0,70);
   hSNR->Draw();
   fitsnr->Draw("lsame");
}

void billr(Int_t compress) {
   //read N histograms from keys
   timer.Start();
   TFile f("bill.root");
   TIter next(f.GetListOfKeys());
   TH1F *h;
   TH1::AddDirectory(kFALSE);
   TKey *key;
   Int_t i=0;
   TH1F *hmean = new TH1F("hmean","hist mean from keys",100,0,1);
   
   while ((key=(TKey*)next())) {
      h = (TH1F*)key->ReadObj();
      hmean->Fill(h->GetMean());
      delete h;
      i++;
   }
   timer.Stop();
   printf("billr%d  : RT=%7.3f s, Cpu=%7.3f s\n",compress,timer.RealTime(),timer.CpuTime());
}

void makePlot(vector<TTree*> sigTree,vector<double> sigWeight,
	      vector<int> ptHatLo, vector<int> ptHatHi,
              std::string var,TCut cut,TH1F* h,bool norm)
{
   TH1F *hRes = (TH1F*)h->Clone();
   hRes->SetName("hRes");
   hRes->Sumw2();
  
   char tmp[300];
   for (unsigned int i=0; i<sigTree.size(); i++)
     {
       // first determine the pthat cut
       
       sprintf(tmp, "ptHat >= %d && ptHat <= %d",ptHatLo[i], ptHatHi[i]);
       TCut ptHatCut = tmp;
       TCut allCut = cut + ptHatCut;
       cout << "Current cut = " << allCut.GetTitle() << endl;
       TH1F *htmp = (TH1F*)h->Clone();
       htmp->SetName("htmp");
       sigTree[i]->Draw(Form("%s>>htmp",var.data()),allCut);
       htmp->Sumw2();
       htmp->Scale(sigWeight[i]);
       cout << "scale = " << sigWeight[i] << endl;
       cout << "After scaling htmp -> entries() " << htmp->GetEntries() << endl;
       cout << "After scaling htmp -> Integral() " << htmp->Integral() << endl;
       cout << "After scaling htmp -> GetMean()  " << htmp->GetMean() << endl;
       cout << "After scaling htmp -> GetRMS()  " << htmp->GetRMS() << endl;
       hRes->Add(htmp);
       delete htmp;
   }
   h->Sumw2();
   h->Add(hRes);
   if(norm)h->Scale(1.0/(double)h->Integral(0,1000));
   cout << "After scaling h-> entries() " << h->GetEntries() << endl;
   cout << "After scaling h-> Integral() " << h->Integral() << endl;
   cout << "After scaling h -> GetMean()  " << h->GetMean() << endl;
   cout << "After scaling h -> GetRMS()  " << h->GetRMS() << endl;
   
   delete hRes;
}

void langaus() {
  // Fill Histogram
  
  TFile *f=new TFile("fitInputs.root");
  TH1F *hSNR = (TH1F*)f->Get("hxsobs");

  // Fitting SNR histo
  printf("Fitting...\n");

  // Setting fit range and start values
  Double_t fr[2];
  Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
  fr[0]=0.3*hSNR->GetMean();
  fr[1]=3.0*hSNR->GetMean();

  pllo[0]=0.5; pllo[1]=5.0; pllo[2]=1.0; pllo[3]=0.4;
  plhi[0]=5.0; plhi[1]=50.0; plhi[2]=1000000.0; plhi[3]=5.0;
  sv[0]=1.8; sv[1]=20.0; sv[2]=50000.0; sv[3]=3.0;

  Double_t chisqr;
  Int_t    ndf;
  TF1 *fitsnr = langaufit(hSNR,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
   
  Double_t SNRPeak, SNRFWHM;
  langaupro(fp,SNRPeak,SNRFWHM);

  printf("Fitting done\nPlotting results...\n");

  // Global style settings
  gStyle->SetOptStat(1111);
  gStyle->SetOptFit(111);
  gStyle->SetLabelSize(0.03,"x");
  gStyle->SetLabelSize(0.03,"y");

  hSNR->GetXaxis()->SetRange(0,70);
  hSNR->Draw();
  fitsnr->Draw("lsame");
}

void DumpRecenterParv2(){
    TString pro = "Pro104";
    int taxi = 8583;
    TFile *fin;
    int nrun = GetTotalRun();
    if(nrun<0) exit("Empty run list file!");

    TFile *fout = new TFile(Form("Recentering.root"),"Recreate");
    TVectorD vecmean;
    TVectorD vecrms;
    vecmean.ResizeTo(ncent*nbbcz*nhar*nsub*nxy);
    vecrms.ResizeTo(ncent*nbbcz*nhar*nsub*nxy);
    for(int irun=0;irun<nrun;irun++){
     cout<<irun<<" of total "<<nrun<<" runs"<<endl;
     int RunNumber=GetRun(irun);
     fin = TFile::Open(Form("Run15pAu200MinBias/output_fvtxwithcntrecenter_%d.root",RunNumber));
     //ofstream fout(Form("Calibration/%s/Recentering_%d.dat",dataset.Data(),GetRun(irun)));
     for(int icent=0;icent<ncent;icent++){
      for(int ibbcz=0;ibbcz<nbbcz;ibbcz++){
       for(int ihar=0;ihar<nhar;ihar++){
        for(int isub=0;isub<nsub;isub++){
         for(int ixy=0;ixy<nxy;ixy++){
             TH1F* q = (TH1F*)fin->Get(Form("q_%d_%d_%d_%d_%d",icent,ibbcz,ihar,isub,ixy));
             float mean = q->GetMean();
             float rms = q->GetRMS();
             vecmean[icent*nbbcz*nhar*nsub*nxy+ibbcz*nhar*nsub*nxy+ihar*nsub*nxy+isub*nxy+ixy] = mean;
             vecrms[icent*nbbcz*nhar*nsub*nxy+ibbcz*nhar*nsub*nxy+ihar*nsub*nxy+isub*nxy+ixy] = rms;
         }
        }
       }
      }
     }
    fout->cd();
    vecmean.Write(Form("mean_%d",RunNumber));
    vecrms.Write(Form("rms_%d",RunNumber));
    fin->Close();
    }
}

void mkChanCorrPlots( void )
{
    Float_t meancopo[kNfls];
    Float_t meancrspo[kNfls];
    Float_t rmscopo[kNfls];
    Float_t rmscrspo[kNfls];

    TFile* fout = TFile::Open( outfn, "recreate" );

    TGraphErrors* gCopo02 = new TGraphErrors;
    TGraphErrors* gCopo13 = new TGraphErrors;
    
    TGraphErrors* gCrspo02 = new TGraphErrors;
    TGraphErrors* gCrspo13 = new TGraphErrors;

    for( UChar_t f = 0; f < kNfls; f++ )
    {
        Float_t mean[kNChCmb];
        Float_t rms[kNChCmb];

        TFile* ifl = TFile::Open( infn[f] );
 
        for( UChar_t c = 0; c < NSnConstants::kNchans; c++ )
        {
            for( UChar_t c2 = c + 1; c2 < NSnConstants::kNchans; c2++ )
            {
                TH1F* h = dynamic_cast<TH1F*>(ifl->Get(Form("hCC%d%d", c, c2)));

                UInt_t idx  = TSnRecoChanOffsets::IndexFor(c2, c);
                mean[idx]   = h->GetMean( );
                rms[idx]    = h->GetRMS( );
            }
        }

        meancopo[f]     = mean[TSnRecoChanOffsets::IndexFor(copoC[f] + 2, 
                                                                   copoC[f])];
        meancrspo[f]    = mean[TSnRecoChanOffsets::IndexFor(crspoC[f] + 2, 
                                                                   crspoC[f])];

        rmscopo[f]      = rms[TSnRecoChanOffsets::IndexFor(copoC[f] + 2, 
                                                                   copoC[f])];
        rmscrspo[f]     = rms[TSnRecoChanOffsets::IndexFor(crspoC[f] + 2, 
                                                                   crspoC[f])];

        Printf( "co po mean: %f rms: %f run: %u", meancopo[f], rmscopo[f], 
                                                                    y[f] );
        Printf( "cross po mean: %f rms: %f run: %u\n", meancrspo[f], 
                                                        rmscrspo[f], y[f] );

        ifl->Close( );
    }

//    Float_t y[kNfls]    = { 1., 1., 1., 1., 1., 1., 1., 1. };
    Float_t ye[kNfls]   = { 0., 0., 0., 0., 0., 0., 0., 0. }; 

    fout->cd( );

    UInt_t nc0 = 0;
    UInt_t nc1 = 0;

    for( UChar_t f = 0; f < kNfls; f++ )
    {
        if( copoC[f] == 0 && crspoC[f] == 1 )
        {
            gCopo02->SetPoint( nc0, meancopo[f], 0. );
            gCrspo13->SetPoint( nc0, meancrspo[f], 0. );

            gCrspo13->SetPointError( nc0, rmscrspo[f], ye[f] );
            gCopo02->SetPointError( nc0++, rmscopo[f], ye[f] );

            gCopo02->SetMarkerColor( kRed );
            gCopo02->SetLineColor( kRed );
            gCopo02->GetHistogram()->SetMinimum( -0.5 );
            gCopo02->GetHistogram()->SetMaximum( 1.5 );

            gCrspo13->SetMarkerColor( kRed );
            gCrspo13->SetLineColor( kRed );
            gCrspo13->GetHistogram()->SetMinimum( -0.5 );
            gCrspo13->GetHistogram()->SetMaximum( 1.5 );
        }            
        if( copoC[f] == 1 && crspoC[f] == 0 )
        {
            gCopo13->SetPoint( nc1, meancopo[f], 1. );
            gCrspo02->SetPoint( nc1, meancrspo[f], 1. );

            gCrspo02->SetPointError( nc1, rmscrspo[f], ye[f] );
            gCopo13->SetPointError( nc1++, rmscopo[f], ye[f] );

            gCopo13->SetMarkerColor( kBlack );
            gCopo13->SetLineColor( kBlack );
            gCopo13->GetHistogram()->SetMinimum( -0.5 );
            gCopo13->GetHistogram()->SetMaximum( 1.5 );

            gCrspo02->SetMarkerColor( kBlack );
            gCrspo02->SetLineColor( kBlack );
            gCrspo02->GetHistogram()->SetMinimum( -0.5 );
            gCrspo02->GetHistogram()->SetMaximum( 1.5 );
        }
    }

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