Jeng/scripts/pvHistory.C

// [USAGE]
// Modify fillNumbers[] to corresponding runs
// .L pvHistory.C+
// pvHistory("skimmed files",true) ==> the second argument means drawing w.r.t LHC fill number instead of runnumber

#if !defined(__CINT__) && !defined(__MAKECINT__)
#include "DataFormats/FWLite/interface/Handle.h"
#include "DataFormats/FWLite/interface/Event.h"
#include "DataFormats/FWLite/interface/ChainEvent.h"
//Headers for the data items
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/BeamSpot/interface/BeamSpot.h"

#endif

#include "TROOT.h"
#include "TH2F.h"
#include "TFile.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TBox.h"
#include <vector>
#include <map>
#include <string>
#include "TStyle.h"
#include "TPaveText.h"
#include "TLegend.h"

using namespace std;
typedef map<unsigned int, unsigned int> RunNumberMap;
typedef map<RunNumberMap,unsigned int> FillNumberMap;

bool LHCFill;
//
// iterative fit of +- 2 sigma Gaussian core on a 1d histogram
//
bool fitHistogram (TH1* slice, TF1& fg, double* pars, double* epars, bool dbg = false) {
  // initialization of parameters
  pars[0] = slice->GetMaximum();
  pars[1] = slice->GetMean();
  pars[2] = slice->GetRMS();
  epars[0] = 0;
  epars[1] = 0;
  epars[2] = 0;
  // some iterations
  bool converged(false);
  for ( unsigned int it=0; it<10; ++it) {
    const double nSigRange = 2.;
    fg.SetParameters(pars);
    double oldMean = pars[1];
    double oldWidth = pars[2];
    if ( it>0 ) {
      fg.SetRange(pars[1]-nSigRange*pars[2],
                  pars[1]+nSigRange*pars[2]);
      //       if ( dbg ) {
      //        cout << "Setting range to " << pars[1]-nSigRange*pars[2]
      //             << " " << pars[1]+nSigRange*pars[2] << endl;
      //       }
    }
    else {
      fg.SetRange(slice->GetXaxis()->GetXmin(),
                  slice->GetXaxis()->GetXmax());
      //       if ( dbg ) {
      //        cout << "Setting range to " << slice->GetXaxis()->GetXmin()
      //             << " " << slice->GetXaxis()->GetXmax() << endl;
      //       }
    }
    if ( dbg ) 
      slice->Fit(&fg,"QLR+");
    else
      slice->Fit(&fg,"QLR");
    fg.GetParameters(pars);
    epars[0] = fg.GetParError(0);
    epars[1] = fg.GetParError(1);
    epars[2] = fg.GetParError(2);
    if ( converged ) {
      cout << "terminating at iteration " << it << " "
           << pars[1] << " " << oldMean << " " << epars[1] << " "
           << pars[2] << " " << oldWidth << " " << epars[2] << endl;
      break;
    }
    else if ( dbg ) {
      cout << "  at iteration " << it << " "
           << pars[1] << " " << oldMean << " " << epars[1] << " "
           << pars[2] << " " << oldWidth << " " << epars[2] << endl;
    }
    // convergence criterion
    if ( it>2 &&
         fabs((pars[1]-oldMean))<max(fabs(0.01*epars[1]),0.0010) &&
         fabs((pars[2]-oldWidth))<max(fabs(0.01*epars[2]),0.0010) ) {
      converged = true;
    }
  }

  return converged;
}

//
// fit Gaussian to slices of the 2d "run" histogram
//
void fitRunHistogram (TH2* histo, TH1*& hmean, TH1*& hwidth) {
  //
  // create histograms for mean and width vs. run
  //
  TDirectory* saveDir = gDirectory;
  gROOT->cd();
  int nbin = histo->GetNbinsX();
  string title = histo->GetName();
  TAxis* xhisto = histo->GetXaxis();
  hmean = new TH1F((title+"Mean").c_str(),title.c_str(),nbin,-0.5,nbin-0.5);
  TAxis* xhmean = hmean->GetXaxis();
  hwidth = new TH1F((title+"Width").c_str(),title.c_str(),nbin,-0.5,nbin-0.5);
  TAxis* xhwidth = hwidth->GetXaxis();
  //
  // fit function
  //
  TF1 fg("fg","gaus");
  //
  // loop over runs
  //
  bool converged;
  double pars[3];
  double epars[3];
  for ( unsigned int ib=1; ib<=nbin; ++ib ) {
    //
    // initialize mean and width to 0 error
    //    
    hmean->SetBinError(ib,0.);
    hwidth->SetBinError(ib,0.);
    //
    // get slice
    //
    TH1* slice = 
      histo->ProjectionY((title+xhisto->GetBinLabel(ib)).c_str(),ib,ib);
    // drop histograms with too few entries
    if ( slice->GetEntries()<100 ) {
      delete slice;
      continue;
    }
    //
    // fit slice
    //
    converged = fitHistogram(slice,fg,pars,epars);
    if ( converged ) {
      // copy parameters and errors
      hmean->SetBinContent(ib,pars[1]);
      hmean->SetBinError(ib,epars[1]);
      hwidth->SetBinContent(ib,pars[2]);
      hwidth->SetBinError(ib,epars[2]);
    }
    else {
      cout << "No convergence for " << histo->GetName() << " / run " 
           << histo->GetXaxis()->GetBinLabel(ib) << endl;
      // refit with debug flag
      // fitHistogram(slice,fg,pars,epars,true);
      hmean->SetBinContent(ib,pars[1]);
      hmean->SetBinError(ib,epars[1]);
      hwidth->SetBinContent(ib,pars[2]);
      hwidth->SetBinError(ib,epars[2]);
    }
    // keep failing fits for debugging
    //     if ( converged )  delete slice;
    xhmean->SetBinLabel(ib,xhisto->GetBinLabel(ib));
    xhwidth->SetBinLabel(ib,xhisto->GetBinLabel(ib));
  }

  saveDir->cd();
}

//
// set x-axis labels to run number
//
void setRunNumberLabels (TH1* histo, RunNumberMap& runNumbers) {
  char label[32];
  TAxis* axis = histo->GetXaxis();
  for ( RunNumberMap::const_iterator i=runNumbers.begin();
        i!=runNumbers.end(); ++i ) {
    if ((*i).first == 124120 || (*i).first == 124275)
      sprintf(label,"%d %s",(*i).first," (2.36 TeV)");
    else
      sprintf(label,"%d",(*i).first);
    axis->SetBinLabel((*i).second+1,label);
  }
}

void setRunNumberLabels (TH1* histo, FillNumberMap& fillNumbers) {
  char label[32];
  TAxis* axis = histo->GetXaxis();
  for ( FillNumberMap::const_iterator i=fillNumbers.begin();
        i!=fillNumbers.end(); ++i ) {
    RunNumberMap tmp((*i).first);
    if (tmp.begin()->first == 124120 || tmp.begin()->first == 124275)
      sprintf(label,"%d %s",(*i).second," (2.36 TeV)");
    else
      sprintf(label,"%d",(*i).second);
    axis->SetBinLabel((tmp.begin()->second)+1,label);
  }
}

//
// plot PV distribution for a single run with indication of BS position
//
void plotRun (char coord, int run) 
{
  char title[64];

  sprintf(title,"pv%cVsRun%d",coord,run);
  TH1* hpv = (TH1*)gROOT->FindObject(title);
  if ( hpv==0 ) return;

  sprintf(title,"pv%cRun%d",coord,run);
  TCanvas* c = new TCanvas(title);

  hpv->Draw();

  sprintf(title,"bs%cVsRunMean",coord);
  TH1* hbs = (TH1*)gROOT->FindObject(title);
  sprintf(title,"%d",run);
  double vbs(0.);
  double ebs(0.);
  for ( unsigned int ib=1; ib<=hbs->GetNbinsX(); ++ib ) {
    const char* label = hbs->GetXaxis()->GetBinLabel(ib);
    if ( strcmp(title,label)==0 ) {
      vbs = hbs->GetBinContent(ib);
      ebs = hbs->GetBinError(ib);
      break;
    }
  }
  TBox* box = new TBox(vbs-ebs,0,vbs+ebs,hpv->GetMaximum());
  box->SetFillStyle(1001);
  box->SetFillColor(7);
  box->Draw();

  hpv->Draw("same");

  sprintf(title,"pv%cRun%d.eps",coord,run);
  c->SaveAs(title);
  delete c;
}
//
// summary plots (PV + BS) vs run
//
void plotPvVsBs () {
  char coord[] = { "xyz" };
  for ( int ic=0; ic<3; ++ic ) {
    string canvasTitle("pvVsBs");
    canvasTitle += coord[ic];
    TCanvas* c = new TCanvas(canvasTitle.c_str(),canvasTitle.c_str(),
                             700,500);
    string bsTitle("bs.VsRunMean");
    bsTitle.replace(2,1,1,coord[ic]);
    TH1* bsVsRun = (TH1*)gROOT->FindObject(bsTitle.c_str());
    string pvTitle("pv.VsRunMean");
    pvTitle.replace(2,1,1,coord[ic]);
    TH1* pvVsRun = (TH1*)gROOT->FindObject(pvTitle.c_str());

    Double_t yMin,yMax;
    if (ic == 2){
      yMax=pvVsRun->GetMaximum()+0.4;
      yMin=pvVsRun->GetMinimum()-0.4;
    }
    else {
      yMax=pvVsRun->GetMaximum()+0.01;
      yMin=pvVsRun->GetMinimum()-0.005;
    }

    pvVsRun->GetYaxis()->SetRangeUser(yMin,yMax);
    pvVsRun->GetXaxis()->LabelsOption("u");
    if (LHCFill)
      pvVsRun->GetXaxis()->SetTitle("LHC Fill");
    else
      pvVsRun->GetXaxis()->SetTitle("Run");
    pvVsRun->GetYaxis()->SetTitle(TString(coord[ic])+" (cm)");
    pvVsRun->Draw();
    bsVsRun->SetFillStyle(1001);
    bsVsRun->SetFillColor(7);
    //     bsVsRun->SetMarkerStyle(21);
    bsVsRun->SetMarkerColor(7);
    bsVsRun->SetLineColor(7);
    bsVsRun->Draw("E2 same");
    pvVsRun->Draw("same");

    TPaveText *ppt= new TPaveText(.42,.7,.82,.86,"NDC");
    ppt->AddText("CMS Preliminary 2009");
    ppt->SetTextAlign(12);
    ppt->SetTextFont(62);
    ppt->SetTextSize(0.045);
    ppt->SetBorderSize(0);
    ppt->SetFillStyle(0);
    ppt->Draw();

    TLegend*l0=new TLegend(0.45,0.6,0.65,0.7);
    l0->SetFillColor(0);
    l0->SetBorderSize(0);
    l0->AddEntry(pvVsRun,"  PrimaryVertex","p");
    l0->AddEntry(bsVsRun,"  Beam Spot","f");
    l0->Draw();

    c->Print(TString(canvasTitle)+".eps");

  }
}

//
// main entry point
//
void pvHistory (const vector<string>& fileNames, bool isFill=false) {
  LHCFill = isFill;

  gROOT->SetStyle("CMS");
  fwlite::ChainEvent ev(fileNames);

  fwlite::Handle< reco::BeamSpot > bsHandle;
  fwlite::Handle< vector<reco::Vertex> > pvHandle;
  //
  // create list of runs with > 100 PVs
  //
  RunNumberMap runNumberCounts;
  int fillNumbers [] = {904,907,907,911,911,911,912,912,916,919,923};
  int lastRun(0);
  for( ev.toBegin(); !ev.atEnd(); ++ev) {
    unsigned int run = ev.id().run();
    pvHandle.getByLabel(ev,"offlinePrimaryVertices");

    const reco::Vertex& pv = (*pvHandle)[0];
    if ( pv.tracksSize()>0 && pv.ndof()>2 && 
         (pv.ndof()+3.)/2./pv.tracksSize()>0.5 ) {
      ++runNumberCounts[run];
    }
  }
  // create list of bin indices
  unsigned int nRun(0);
  RunNumberMap runNumberIndices;
  FillNumberMap fillNumberIndices;
  int nFill=0;
  for ( RunNumberMap::iterator i=runNumberCounts.begin();
        i!=runNumberCounts.end(); ++i,nRun++,nFill++ ) {
    if ( (*i).second>100 ) {
      runNumberIndices[(*i).first] = nRun;
      RunNumberMap tmp;
      tmp[(*i).first] = nRun;
      fillNumberIndices[tmp] = fillNumbers[nFill];
    }
  }
  for ( FillNumberMap::iterator i=fillNumberIndices.begin();
        i!=fillNumberIndices.end();++i){
    RunNumberMap tmp((*i).first);
    cout << "Run = " << tmp.begin()->first << "; Fill = " << (*i).second << endl;
  }
  cout << "Found " << nRun << " run numbers" << endl;
  //
  // create histograms
  // 
  TDirectory* saveDir = gDirectory;
  gROOT->cd();
  TH2* pvxVsRun = new TH2F("pvxVsRun","pvxVsRun",nRun,-0.5,nRun-0.5,80,-0.1,0.3);
  TH2* pvyVsRun = new TH2F("pvyVsRun","pvyVsRun",nRun,-0.5,nRun-0.5,80,-0.2,0.2);
  TH2* pvzVsRun = new TH2F("pvzVsRun","pvzVsRun",nRun,-0.5,nRun-0.5,80,-20,20);
  TH2* bsxVsRun = new TH2F("bsxVsRun","bsxVsRun",nRun,-0.5,nRun-0.5,80,0,0.4);
  TH2* bsyVsRun = new TH2F("bsyVsRun","bsyVsRun",nRun,-0.5,nRun-0.5,80,0,0.4);
  TH2* bszVsRun = new TH2F("bszVsRun","bszVsRun",nRun,-0.5,nRun-0.5,80,-20,20);
  TH1* bsxVsRunMean = new TH1F("bsxVsRunMean","bsxVsRunMean",nRun,-0.5,nRun-0.5);
  TH1* bsyVsRunMean = new TH1F("bsyVsRunMean","bsyVsRunMean",nRun,-0.5,nRun-0.5);
  TH1* bszVsRunMean = new TH1F("bszVsRunMean","bszVsRunMean",nRun,-0.5,nRun-0.5);
  saveDir->cd();

  //
  // fill histograms
  //
  lastRun = -1;
  int idx(-1);
  //int nev(0);
  for( ev.toBegin(); !ev.atEnd(); ++ev) {
    int run = ev.id().run();
    //     if ( (++nev)%1000 == 0 )  cout << "processing record " << nev << " ( run "
    //                            << ev.id().run() << " " << " evt "
    //                            << ev.id().event() << " )" << endl;
    //
    // update bin index if run changed
    //
    if ( run!=lastRun ) {
      RunNumberMap::const_iterator irun = runNumberIndices.find(run);
      if ( irun == runNumberIndices.end() ) {
        idx = -1;
      }
      else {
        idx = (*irun).second;
      }
      lastRun = run;
      //       cout << "Run " << run << " idx " << idx << endl;
    }
    if ( idx<0 )  continue;

    //
    // retrieve beam spot and PV
    //
    bsHandle.getByLabel(ev,"offlineBeamSpot");
    pvHandle.getByLabel(ev,"offlinePrimaryVertices");
    //
    // fill PV histograms (after quality cut)
    //
    const reco::Vertex& pv = (*pvHandle)[0];
    if ( pv.tracksSize()>0 && pv.ndof()>2 && 
         (pv.ndof()+3.)/2./pv.tracksSize()>0.5 ) {
      pvxVsRun->Fill(idx,pv.x());
      pvyVsRun->Fill(idx,pv.y());
      pvzVsRun->Fill(idx,pv.z());
    }
    //
    // fill BS histograms
    //
    bsxVsRun->Fill(idx,bsHandle->x0());
    bsyVsRun->Fill(idx,bsHandle->y0());
    bszVsRun->Fill(idx,bsHandle->z0());
    if ( idx>=0 && idx<bsxVsRun->GetNbinsX() ) {
      // assume BS is constant / run: fill current value into 1d histogram
      bsxVsRunMean->SetBinContent(idx+1,bsHandle->x0());
      bsxVsRunMean->SetBinError(idx+1,bsHandle->x0Error());
      bsyVsRunMean->SetBinContent(idx+1,bsHandle->y0());
      bsyVsRunMean->SetBinError(idx+1,bsHandle->y0Error());
      bszVsRunMean->SetBinContent(idx+1,bsHandle->z0());
      bszVsRunMean->SetBinError(idx+1,bsHandle->z0Error());
    }

  }
  saveDir = gDirectory;
  //
  // number of entries / run
  //
  gROOT->cd();
  TH1* npvVsRun = pvzVsRun->ProjectionX("npvVsRun");
  npvVsRun->SetTitle("npvVsRun");
  TH1* nbsVsRun = bszVsRun->ProjectionX("nbsVsRun");
  npvVsRun->SetTitle("nbsVsRun");
  saveDir->cd();
  //
  // set run numbers on x-axis labels
  //
  if (! LHCFill) {
    setRunNumberLabels(npvVsRun,runNumberIndices);
    setRunNumberLabels(pvxVsRun,runNumberIndices);
    setRunNumberLabels(pvyVsRun,runNumberIndices);
    setRunNumberLabels(pvzVsRun,runNumberIndices);
    setRunNumberLabels(nbsVsRun,runNumberIndices);
    setRunNumberLabels(bsxVsRun,runNumberIndices);
    setRunNumberLabels(bsyVsRun,runNumberIndices);
    setRunNumberLabels(bszVsRun,runNumberIndices);
    setRunNumberLabels(bsxVsRunMean,runNumberIndices);
    setRunNumberLabels(bsyVsRunMean,runNumberIndices);
    setRunNumberLabels(bszVsRunMean,runNumberIndices);
  }
  else{
    setRunNumberLabels(npvVsRun,fillNumberIndices);
    setRunNumberLabels(pvxVsRun,fillNumberIndices);
    setRunNumberLabels(pvyVsRun,fillNumberIndices);
    setRunNumberLabels(pvzVsRun,fillNumberIndices);
    setRunNumberLabels(nbsVsRun,fillNumberIndices);
    setRunNumberLabels(bsxVsRun,fillNumberIndices);
    setRunNumberLabels(bsyVsRun,fillNumberIndices);
    setRunNumberLabels(bszVsRun,fillNumberIndices);
    setRunNumberLabels(bsxVsRunMean,fillNumberIndices);
    setRunNumberLabels(bsyVsRunMean,fillNumberIndices);
    setRunNumberLabels(bszVsRunMean,fillNumberIndices);
  }
  //
  // fit histograms
  //
  TH1* pvxMean; TH1* pvxWidth;
  fitRunHistogram(pvxVsRun,pvxMean,pvxWidth);
  TH1* pvyMean; TH1* pvyWidth;
  fitRunHistogram(pvyVsRun,pvyMean,pvyWidth);
  TH1* pvzMean; TH1* pvzWidth;
  fitRunHistogram(pvzVsRun,pvzMean,pvzWidth);

  //   cout << pvxMean << endl;
  //   cout << pvxMean->GetNbinsX() << " " << pvxMean->GetXaxis() << endl;
  //
  // summary table
  // 
  for ( unsigned int i=1; i<=pvxMean->GetNbinsX(); ++i ) {
    if ( pvxMean->GetBinError(i)>1.e-10 || 
         pvxMean->GetBinError(i)>1.e-10 || 
         pvxMean->GetBinError(i)>1.e-10 ) {
      printf(" %-6s :  %8.4f +- %6.4f  %6.4f +- %6.4f   %8.4f +- %6.4f  %6.4f +- %6.4f   %9.4f +- %7.4f  %6.4f +- %7.4f\n",
             pvxMean->GetXaxis()->GetBinLabel(i),
             pvxMean->GetBinContent(i),pvxMean->GetBinError(i),
             pvxWidth->GetBinContent(i),pvxWidth->GetBinError(i),
             pvyMean->GetBinContent(i),pvyMean->GetBinError(i),
             pvyWidth->GetBinContent(i),pvyWidth->GetBinError(i),
             pvzMean->GetBinContent(i),pvzMean->GetBinError(i),
             pvzWidth->GetBinContent(i),pvzWidth->GetBinError(i));
    }
  }

  //
  // save histograms
  //
  TFile* fo = new TFile("./skim_BSC0_AND_40OR41_MinBias_Dec19th_v1_plots.root","recreate");
  pvxVsRun->Write();
  pvyVsRun->Write();
  pvzVsRun->Write();
  pvxMean->Write();
  pvxWidth->Write();
  pvyMean->Write();
  pvyWidth->Write();
  pvzMean->Write();
  pvzWidth->Write();
  fo->Write();
  delete fo;

  plotPvVsBs();
  //
  // plot inidividual histograms / run
  //
  //   for ( RunNumberMap::iterator i=runNumberIndices.begin();
  //    i!=runNumberIndices.end(); ++i ) {
  //     plotRun('x',(*i).first);
  //     plotRun('y',(*i).first);
  //     plotRun('z',(*i).first);
  //   }  

}
//
// call with single file name
//
void pvHistory (const char* fileName) {
  vector<string> fileNames(1,fileName);
  pvHistory(fileNames,false);
}

void pvHistory (const char* fileName,bool isFill) {
  vector<string> fileNames(1,fileName);
  pvHistory(fileNames,isFill);
}

Revision:	1.4
Committed:	Thu Feb 18 00:38:45 2010 UTC (15 years, 2 months ago) by jengbou
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.3:	+4 -1 lines
Log Message:	Fix x label
#	User	Rev	Content
1	jengbou	1.3	// [USAGE]
2			// Modify fillNumbers[] to corresponding runs
3			// .L pvHistory.C+
4			// pvHistory("skimmed files",true) ==> the second argument means drawing w.r.t LHC fill number instead of runnumber
5
6	jengbou	1.1	#if !defined(__CINT__) && !defined(__MAKECINT__)
7			#include "DataFormats/FWLite/interface/Handle.h"
8			#include "DataFormats/FWLite/interface/Event.h"
9			#include "DataFormats/FWLite/interface/ChainEvent.h"
10			//Headers for the data items
11			#include "DataFormats/VertexReco/interface/Vertex.h"
12			#include "DataFormats/BeamSpot/interface/BeamSpot.h"
13
14			#endif
15
16			#include "TROOT.h"
17			#include "TH2F.h"
18			#include "TFile.h"
19			#include "TF1.h"
20			#include "TCanvas.h"
21			#include "TMath.h"
22			#include "TBox.h"
23			#include <vector>
24			#include <map>
25			#include <string>
26	jengbou	1.2	#include "TStyle.h"
27			#include "TPaveText.h"
28			#include "TLegend.h"
29	jengbou	1.1
30			using namespace std;
31			typedef map<unsigned int, unsigned int> RunNumberMap;
32	jengbou	1.3	typedef map<RunNumberMap,unsigned int> FillNumberMap;
33	jengbou	1.1
34	jengbou	1.3	bool LHCFill;
35	jengbou	1.1	//
36			// iterative fit of +- 2 sigma Gaussian core on a 1d histogram
37			//
38			bool fitHistogram (TH1* slice, TF1& fg, double* pars, double* epars, bool dbg = false) {
39			// initialization of parameters
40			pars[0] = slice->GetMaximum();
41			pars[1] = slice->GetMean();
42			pars[2] = slice->GetRMS();
43			epars[0] = 0;
44			epars[1] = 0;
45			epars[2] = 0;
46			// some iterations
47			bool converged(false);
48			for ( unsigned int it=0; it<10; ++it) {
49			const double nSigRange = 2.;
50			fg.SetParameters(pars);
51			double oldMean = pars[1];
52			double oldWidth = pars[2];
53			if ( it>0 ) {
54			fg.SetRange(pars[1]-nSigRange*pars[2],
55			pars[1]+nSigRange*pars[2]);
56	jengbou	1.2	// if ( dbg ) {
57			// cout << "Setting range to " << pars[1]-nSigRange*pars[2]
58			// << " " << pars[1]+nSigRange*pars[2] << endl;
59			// }
60	jengbou	1.1	}
61			else {
62			fg.SetRange(slice->GetXaxis()->GetXmin(),
63			slice->GetXaxis()->GetXmax());
64	jengbou	1.2	// if ( dbg ) {
65			// cout << "Setting range to " << slice->GetXaxis()->GetXmin()
66			// << " " << slice->GetXaxis()->GetXmax() << endl;
67			// }
68	jengbou	1.1	}
69			if ( dbg )
70			slice->Fit(&fg,"QLR+");
71			else
72			slice->Fit(&fg,"QLR");
73			fg.GetParameters(pars);
74			epars[0] = fg.GetParError(0);
75			epars[1] = fg.GetParError(1);
76			epars[2] = fg.GetParError(2);
77			if ( converged ) {
78			cout << "terminating at iteration " << it << " "
79			<< pars[1] << " " << oldMean << " " << epars[1] << " "
80			<< pars[2] << " " << oldWidth << " " << epars[2] << endl;
81			break;
82			}
83			else if ( dbg ) {
84			cout << " at iteration " << it << " "
85			<< pars[1] << " " << oldMean << " " << epars[1] << " "
86			<< pars[2] << " " << oldWidth << " " << epars[2] << endl;
87			}
88			// convergence criterion
89			if ( it>2 &&
90			fabs((pars[1]-oldMean))<max(fabs(0.01*epars[1]),0.0010) &&
91			fabs((pars[2]-oldWidth))<max(fabs(0.01*epars[2]),0.0010) ) {
92			converged = true;
93			}
94			}
95
96			return converged;
97			}
98
99			//
100			// fit Gaussian to slices of the 2d "run" histogram
101			//
102			void fitRunHistogram (TH2* histo, TH1& hmean, TH1& hwidth) {
103			//
104			// create histograms for mean and width vs. run
105			//
106			TDirectory* saveDir = gDirectory;
107			gROOT->cd();
108			int nbin = histo->GetNbinsX();
109			string title = histo->GetName();
110			TAxis* xhisto = histo->GetXaxis();
111			hmean = new TH1F((title+"Mean").c_str(),title.c_str(),nbin,-0.5,nbin-0.5);
112			TAxis* xhmean = hmean->GetXaxis();
113			hwidth = new TH1F((title+"Width").c_str(),title.c_str(),nbin,-0.5,nbin-0.5);
114			TAxis* xhwidth = hwidth->GetXaxis();
115			//
116			// fit function
117			//
118			TF1 fg("fg","gaus");
119			//
120			// loop over runs
121			//
122			bool converged;
123			double pars[3];
124			double epars[3];
125			for ( unsigned int ib=1; ib<=nbin; ++ib ) {
126			//
127			// initialize mean and width to 0 error
128			//
129			hmean->SetBinError(ib,0.);
130			hwidth->SetBinError(ib,0.);
131			//
132			// get slice
133			//
134			TH1* slice =
135			histo->ProjectionY((title+xhisto->GetBinLabel(ib)).c_str(),ib,ib);
136			// drop histograms with too few entries
137			if ( slice->GetEntries()<100 ) {
138			delete slice;
139			continue;
140			}
141			//
142			// fit slice
143			//
144			converged = fitHistogram(slice,fg,pars,epars);
145			if ( converged ) {
146			// copy parameters and errors
147			hmean->SetBinContent(ib,pars[1]);
148			hmean->SetBinError(ib,epars[1]);
149			hwidth->SetBinContent(ib,pars[2]);
150			hwidth->SetBinError(ib,epars[2]);
151			}
152			else {
153			cout << "No convergence for " << histo->GetName() << " / run "
154			<< histo->GetXaxis()->GetBinLabel(ib) << endl;
155			// refit with debug flag
156			// fitHistogram(slice,fg,pars,epars,true);
157			hmean->SetBinContent(ib,pars[1]);
158			hmean->SetBinError(ib,epars[1]);
159			hwidth->SetBinContent(ib,pars[2]);
160			hwidth->SetBinError(ib,epars[2]);
161			}
162			// keep failing fits for debugging
163	jengbou	1.2	// if ( converged ) delete slice;
164	jengbou	1.1	xhmean->SetBinLabel(ib,xhisto->GetBinLabel(ib));
165			xhwidth->SetBinLabel(ib,xhisto->GetBinLabel(ib));
166			}
167
168			saveDir->cd();
169			}
170
171			//
172			// set x-axis labels to run number
173			//
174			void setRunNumberLabels (TH1* histo, RunNumberMap& runNumbers) {
175			char label[32];
176			TAxis* axis = histo->GetXaxis();
177			for ( RunNumberMap::const_iterator i=runNumbers.begin();
178			i!=runNumbers.end(); ++i ) {
179	jengbou	1.2	if ((i).first == 124120 \|\| (i).first == 124275)
180			sprintf(label,"%d %s",(*i).first," (2.36 TeV)");
181			else
182			sprintf(label,"%d",(*i).first);
183	jengbou	1.1	axis->SetBinLabel((*i).second+1,label);
184			}
185			}
186
187	jengbou	1.3	void setRunNumberLabels (TH1* histo, FillNumberMap& fillNumbers) {
188			char label[32];
189			TAxis* axis = histo->GetXaxis();
190			for ( FillNumberMap::const_iterator i=fillNumbers.begin();
191			i!=fillNumbers.end(); ++i ) {
192			RunNumberMap tmp((*i).first);
193			if (tmp.begin()->first == 124120 \|\| tmp.begin()->first == 124275)
194			sprintf(label,"%d %s",(*i).second," (2.36 TeV)");
195			else
196			sprintf(label,"%d",(*i).second);
197			axis->SetBinLabel((tmp.begin()->second)+1,label);
198			}
199			}
200	jengbou	1.1
201			//
202			// plot PV distribution for a single run with indication of BS position
203			//
204			void plotRun (char coord, int run)
205			{
206			char title[64];
207
208			sprintf(title,"pv%cVsRun%d",coord,run);
209			TH1* hpv = (TH1*)gROOT->FindObject(title);
210			if ( hpv==0 ) return;
211
212			sprintf(title,"pv%cRun%d",coord,run);
213			TCanvas* c = new TCanvas(title);
214
215			hpv->Draw();
216
217			sprintf(title,"bs%cVsRunMean",coord);
218			TH1* hbs = (TH1*)gROOT->FindObject(title);
219			sprintf(title,"%d",run);
220			double vbs(0.);
221			double ebs(0.);
222			for ( unsigned int ib=1; ib<=hbs->GetNbinsX(); ++ib ) {
223			const char* label = hbs->GetXaxis()->GetBinLabel(ib);
224			if ( strcmp(title,label)==0 ) {
225			vbs = hbs->GetBinContent(ib);
226			ebs = hbs->GetBinError(ib);
227			break;
228			}
229			}
230			TBox* box = new TBox(vbs-ebs,0,vbs+ebs,hpv->GetMaximum());
231			box->SetFillStyle(1001);
232			box->SetFillColor(7);
233			box->Draw();
234
235			hpv->Draw("same");
236
237	jengbou	1.3	sprintf(title,"pv%cRun%d.eps",coord,run);
238	jengbou	1.1	c->SaveAs(title);
239			delete c;
240			}
241			//
242			// summary plots (PV + BS) vs run
243			//
244			void plotPvVsBs () {
245			char coord[] = { "xyz" };
246			for ( int ic=0; ic<3; ++ic ) {
247			string canvasTitle("pvVsBs");
248			canvasTitle += coord[ic];
249			TCanvas* c = new TCanvas(canvasTitle.c_str(),canvasTitle.c_str(),
250			700,500);
251			string bsTitle("bs.VsRunMean");
252			bsTitle.replace(2,1,1,coord[ic]);
253			TH1* bsVsRun = (TH1*)gROOT->FindObject(bsTitle.c_str());
254			string pvTitle("pv.VsRunMean");
255			pvTitle.replace(2,1,1,coord[ic]);
256			TH1* pvVsRun = (TH1*)gROOT->FindObject(pvTitle.c_str());
257	jengbou	1.2
258			Double_t yMin,yMax;
259			if (ic == 2){
260			yMax=pvVsRun->GetMaximum()+0.4;
261			yMin=pvVsRun->GetMinimum()-0.4;
262			}
263			else {
264			yMax=pvVsRun->GetMaximum()+0.01;
265			yMin=pvVsRun->GetMinimum()-0.005;
266			}
267
268			pvVsRun->GetYaxis()->SetRangeUser(yMin,yMax);
269			pvVsRun->GetXaxis()->LabelsOption("u");
270	jengbou	1.4	if (LHCFill)
271			pvVsRun->GetXaxis()->SetTitle("LHC Fill");
272			else
273			pvVsRun->GetXaxis()->SetTitle("Run");
274	jengbou	1.2	pvVsRun->GetYaxis()->SetTitle(TString(coord[ic])+" (cm)");
275	jengbou	1.1	pvVsRun->Draw();
276			bsVsRun->SetFillStyle(1001);
277			bsVsRun->SetFillColor(7);
278	jengbou	1.2	// bsVsRun->SetMarkerStyle(21);
279			bsVsRun->SetMarkerColor(7);
280			bsVsRun->SetLineColor(7);
281	jengbou	1.1	bsVsRun->Draw("E2 same");
282			pvVsRun->Draw("same");
283	jengbou	1.2
284			TPaveText *ppt= new TPaveText(.42,.7,.82,.86,"NDC");
285			ppt->AddText("CMS Preliminary 2009");
286			ppt->SetTextAlign(12);
287			ppt->SetTextFont(62);
288			ppt->SetTextSize(0.045);
289			ppt->SetBorderSize(0);
290			ppt->SetFillStyle(0);
291			ppt->Draw();
292
293			TLegend*l0=new TLegend(0.45,0.6,0.65,0.7);
294			l0->SetFillColor(0);
295			l0->SetBorderSize(0);
296			l0->AddEntry(pvVsRun," PrimaryVertex","p");
297			l0->AddEntry(bsVsRun," Beam Spot","f");
298			l0->Draw();
299
300	jengbou	1.3	c->Print(TString(canvasTitle)+".eps");
301	jengbou	1.2
302	jengbou	1.1	}
303			}
304
305			//
306			// main entry point
307			//
308	jengbou	1.3	void pvHistory (const vector<string>& fileNames, bool isFill=false) {
309			LHCFill = isFill;
310	jengbou	1.1
311	jengbou	1.2	gROOT->SetStyle("CMS");
312	jengbou	1.1	fwlite::ChainEvent ev(fileNames);
313
314			fwlite::Handle< reco::BeamSpot > bsHandle;
315			fwlite::Handle< vector<reco::Vertex> > pvHandle;
316			//
317			// create list of runs with > 100 PVs
318			//
319			RunNumberMap runNumberCounts;
320	jengbou	1.3	int fillNumbers [] = {904,907,907,911,911,911,912,912,916,919,923};
321	jengbou	1.1	int lastRun(0);
322			for( ev.toBegin(); !ev.atEnd(); ++ev) {
323			unsigned int run = ev.id().run();
324			pvHandle.getByLabel(ev,"offlinePrimaryVertices");
325
326			const reco::Vertex& pv = (*pvHandle)[0];
327			if ( pv.tracksSize()>0 && pv.ndof()>2 &&
328			(pv.ndof()+3.)/2./pv.tracksSize()>0.5 ) {
329			++runNumberCounts[run];
330			}
331			}
332			// create list of bin indices
333			unsigned int nRun(0);
334			RunNumberMap runNumberIndices;
335	jengbou	1.3	FillNumberMap fillNumberIndices;
336			int nFill=0;
337	jengbou	1.1	for ( RunNumberMap::iterator i=runNumberCounts.begin();
338	jengbou	1.3	i!=runNumberCounts.end(); ++i,nRun++,nFill++ ) {
339			if ( (*i).second>100 ) {
340			runNumberIndices[(*i).first] = nRun;
341			RunNumberMap tmp;
342			tmp[(*i).first] = nRun;
343			fillNumberIndices[tmp] = fillNumbers[nFill];
344			}
345			}
346			for ( FillNumberMap::iterator i=fillNumberIndices.begin();
347			i!=fillNumberIndices.end();++i){
348			RunNumberMap tmp((*i).first);
349			cout << "Run = " << tmp.begin()->first << "; Fill = " << (*i).second << endl;
350	jengbou	1.1	}
351			cout << "Found " << nRun << " run numbers" << endl;
352			//
353			// create histograms
354			//
355			TDirectory* saveDir = gDirectory;
356			gROOT->cd();
357			TH2* pvxVsRun = new TH2F("pvxVsRun","pvxVsRun",nRun,-0.5,nRun-0.5,80,-0.1,0.3);
358			TH2* pvyVsRun = new TH2F("pvyVsRun","pvyVsRun",nRun,-0.5,nRun-0.5,80,-0.2,0.2);
359			TH2* pvzVsRun = new TH2F("pvzVsRun","pvzVsRun",nRun,-0.5,nRun-0.5,80,-20,20);
360			TH2* bsxVsRun = new TH2F("bsxVsRun","bsxVsRun",nRun,-0.5,nRun-0.5,80,0,0.4);
361			TH2* bsyVsRun = new TH2F("bsyVsRun","bsyVsRun",nRun,-0.5,nRun-0.5,80,0,0.4);
362			TH2* bszVsRun = new TH2F("bszVsRun","bszVsRun",nRun,-0.5,nRun-0.5,80,-20,20);
363			TH1* bsxVsRunMean = new TH1F("bsxVsRunMean","bsxVsRunMean",nRun,-0.5,nRun-0.5);
364			TH1* bsyVsRunMean = new TH1F("bsyVsRunMean","bsyVsRunMean",nRun,-0.5,nRun-0.5);
365			TH1* bszVsRunMean = new TH1F("bszVsRunMean","bszVsRunMean",nRun,-0.5,nRun-0.5);
366			saveDir->cd();
367
368			//
369			// fill histograms
370			//
371			lastRun = -1;
372			int idx(-1);
373	jengbou	1.3	//int nev(0);
374	jengbou	1.1	for( ev.toBegin(); !ev.atEnd(); ++ev) {
375			int run = ev.id().run();
376	jengbou	1.2	// if ( (++nev)%1000 == 0 ) cout << "processing record " << nev << " ( run "
377			// << ev.id().run() << " " << " evt "
378			// << ev.id().event() << " )" << endl;
379	jengbou	1.1	//
380			// update bin index if run changed
381			//
382			if ( run!=lastRun ) {
383			RunNumberMap::const_iterator irun = runNumberIndices.find(run);
384			if ( irun == runNumberIndices.end() ) {
385			idx = -1;
386			}
387			else {
388			idx = (*irun).second;
389			}
390			lastRun = run;
391	jengbou	1.2	// cout << "Run " << run << " idx " << idx << endl;
392	jengbou	1.1	}
393			if ( idx<0 ) continue;
394
395			//
396			// retrieve beam spot and PV
397			//
398			bsHandle.getByLabel(ev,"offlineBeamSpot");
399			pvHandle.getByLabel(ev,"offlinePrimaryVertices");
400			//
401			// fill PV histograms (after quality cut)
402			//
403			const reco::Vertex& pv = (*pvHandle)[0];
404			if ( pv.tracksSize()>0 && pv.ndof()>2 &&
405			(pv.ndof()+3.)/2./pv.tracksSize()>0.5 ) {
406			pvxVsRun->Fill(idx,pv.x());
407			pvyVsRun->Fill(idx,pv.y());
408			pvzVsRun->Fill(idx,pv.z());
409			}
410			//
411			// fill BS histograms
412			//
413			bsxVsRun->Fill(idx,bsHandle->x0());
414			bsyVsRun->Fill(idx,bsHandle->y0());
415			bszVsRun->Fill(idx,bsHandle->z0());
416			if ( idx>=0 && idx<bsxVsRun->GetNbinsX() ) {
417			// assume BS is constant / run: fill current value into 1d histogram
418			bsxVsRunMean->SetBinContent(idx+1,bsHandle->x0());
419			bsxVsRunMean->SetBinError(idx+1,bsHandle->x0Error());
420			bsyVsRunMean->SetBinContent(idx+1,bsHandle->y0());
421			bsyVsRunMean->SetBinError(idx+1,bsHandle->y0Error());
422			bszVsRunMean->SetBinContent(idx+1,bsHandle->z0());
423			bszVsRunMean->SetBinError(idx+1,bsHandle->z0Error());
424			}
425
426			}
427			saveDir = gDirectory;
428			//
429			// number of entries / run
430			//
431			gROOT->cd();
432			TH1* npvVsRun = pvzVsRun->ProjectionX("npvVsRun");
433			npvVsRun->SetTitle("npvVsRun");
434			TH1* nbsVsRun = bszVsRun->ProjectionX("nbsVsRun");
435			npvVsRun->SetTitle("nbsVsRun");
436			saveDir->cd();
437			//
438			// set run numbers on x-axis labels
439			//
440	jengbou	1.3	if (! LHCFill) {
441			setRunNumberLabels(npvVsRun,runNumberIndices);
442			setRunNumberLabels(pvxVsRun,runNumberIndices);
443			setRunNumberLabels(pvyVsRun,runNumberIndices);
444			setRunNumberLabels(pvzVsRun,runNumberIndices);
445			setRunNumberLabels(nbsVsRun,runNumberIndices);
446			setRunNumberLabels(bsxVsRun,runNumberIndices);
447			setRunNumberLabels(bsyVsRun,runNumberIndices);
448			setRunNumberLabels(bszVsRun,runNumberIndices);
449			setRunNumberLabels(bsxVsRunMean,runNumberIndices);
450			setRunNumberLabels(bsyVsRunMean,runNumberIndices);
451			setRunNumberLabels(bszVsRunMean,runNumberIndices);
452			}
453			else{
454			setRunNumberLabels(npvVsRun,fillNumberIndices);
455			setRunNumberLabels(pvxVsRun,fillNumberIndices);
456			setRunNumberLabels(pvyVsRun,fillNumberIndices);
457			setRunNumberLabels(pvzVsRun,fillNumberIndices);
458			setRunNumberLabels(nbsVsRun,fillNumberIndices);
459			setRunNumberLabels(bsxVsRun,fillNumberIndices);
460			setRunNumberLabels(bsyVsRun,fillNumberIndices);
461			setRunNumberLabels(bszVsRun,fillNumberIndices);
462			setRunNumberLabels(bsxVsRunMean,fillNumberIndices);
463			setRunNumberLabels(bsyVsRunMean,fillNumberIndices);
464			setRunNumberLabels(bszVsRunMean,fillNumberIndices);
465			}
466	jengbou	1.1	//
467			// fit histograms
468			//
469			TH1* pvxMean; TH1* pvxWidth;
470			fitRunHistogram(pvxVsRun,pvxMean,pvxWidth);
471			TH1* pvyMean; TH1* pvyWidth;
472			fitRunHistogram(pvyVsRun,pvyMean,pvyWidth);
473			TH1* pvzMean; TH1* pvzWidth;
474			fitRunHistogram(pvzVsRun,pvzMean,pvzWidth);
475
476	jengbou	1.2	// cout << pvxMean << endl;
477			// cout << pvxMean->GetNbinsX() << " " << pvxMean->GetXaxis() << endl;
478	jengbou	1.1	//
479			// summary table
480			//
481			for ( unsigned int i=1; i<=pvxMean->GetNbinsX(); ++i ) {
482			if ( pvxMean->GetBinError(i)>1.e-10 \|\|
483			pvxMean->GetBinError(i)>1.e-10 \|\|
484			pvxMean->GetBinError(i)>1.e-10 ) {
485			printf(" %-6s : %8.4f +- %6.4f %6.4f +- %6.4f %8.4f +- %6.4f %6.4f +- %6.4f %9.4f +- %7.4f %6.4f +- %7.4f\n",
486			pvxMean->GetXaxis()->GetBinLabel(i),
487			pvxMean->GetBinContent(i),pvxMean->GetBinError(i),
488			pvxWidth->GetBinContent(i),pvxWidth->GetBinError(i),
489			pvyMean->GetBinContent(i),pvyMean->GetBinError(i),
490			pvyWidth->GetBinContent(i),pvyWidth->GetBinError(i),
491			pvzMean->GetBinContent(i),pvzMean->GetBinError(i),
492			pvzWidth->GetBinContent(i),pvzWidth->GetBinError(i));
493			}
494			}
495
496			//
497			// save histograms
498			//
499	jengbou	1.3	TFile* fo = new TFile("./skim_BSC0_AND_40OR41_MinBias_Dec19th_v1_plots.root","recreate");
500	jengbou	1.1	pvxVsRun->Write();
501			pvyVsRun->Write();
502			pvzVsRun->Write();
503			pvxMean->Write();
504			pvxWidth->Write();
505			pvyMean->Write();
506			pvyWidth->Write();
507			pvzMean->Write();
508			pvzWidth->Write();
509			fo->Write();
510			delete fo;
511
512			plotPvVsBs();
513			//
514			// plot inidividual histograms / run
515			//
516	jengbou	1.2	// for ( RunNumberMap::iterator i=runNumberIndices.begin();
517			// i!=runNumberIndices.end(); ++i ) {
518			// plotRun('x',(*i).first);
519			// plotRun('y',(*i).first);
520			// plotRun('z',(*i).first);
521			// }
522	jengbou	1.1
523			}
524			//
525			// call with single file name
526			//
527			void pvHistory (const char* fileName) {
528			vector<string> fileNames(1,fileName);
529	jengbou	1.3	pvHistory(fileNames,false);
530			}
531
532			void pvHistory (const char* fileName,bool isFill) {
533			vector<string> fileNames(1,fileName);
534			pvHistory(fileNames,isFill);
535	jengbou	1.1	}
536	jengbou	1.2