C++ TTree::Project方法代码示例

void drawHFTime(char file[100]="/uscms_data/d1/sunanda/CMSSW_1_3_3/src/SimG4CMS/Calo/test/vcal5x5.root", int mode=0, float xmax=1000) {

  setTDRStyle();
  TFile *File1 = new TFile(file);
  TTree *em = File1->Get("h3");
  TTree *had = File1->Get("h8");
  TH1F *h1 = new TH1F("EM","Time (EM)",500,0,xmax);
  TH1F *h2 = new TH1F("Had","Time (Had)",500,0,xmax);
  if (mode <=0) {
    em->Project("EM","It");
    had->Project("Had","It");
  } else {
    em->Project("EM","T");
    had->Project("Had","T");
  }
  h1->GetXaxis()->SetTitle("Time (EM)");
  h2->GetXaxis()->SetTitle("Time (Had)");

  TCanvas *em1 = new TCanvas("em1","",800,500);
  gPad->SetLogy(1);
  h1->Draw();

  TCanvas *had1 = new TCanvas("had1","",800,500);
  gPad->SetLogy(1);
  h2->Draw();
}

void jetPlot(std::string runNumber = "276282") {

  std::string fileName = "/data/dasu/l1tCaloSummary-" + runNumber + ".root";

  TFile file1(fileName.c_str());

  TTree *tree = (TTree*) file1.Get("Events");

  TH1F *cJetPT = new TH1F("CentralJetPT", "PT Distribution", 300, 0, 300);
  tree->Project("CentralJetPT", "l1extraL1JetParticles_uct2016EmulatorDigis_Central_L1TCaloSummaryTest.obj.m_state.p4Polar_.fCoordinates.fPt");
  cJetPT->SetLineColor(kBlue);

  TH1F *fJetPT = new TH1F("ForwardJetPT", "PT Distribution", 300, 0, 300);
  tree->Project("ForwardJetPT", "l1extraL1JetParticles_uct2016EmulatorDigis_Forward_L1TCaloSummaryTest.obj.m_state.p4Polar_.fCoordinates.fPt");
  fJetPT->SetLineColor(kRed);


  TCanvas *canvas = new TCanvas();
  gPad->SetLogy(1);
  gStyle->SetOptStat(0);
  cJetPT->Draw("HIST,E1");
  fJetPT->Draw("SAME,HIST,E1");

  TLegend *lg=new TLegend(0.55,0.55,0.85,0.85);
  lg->SetFillColor(kWhite);
  lg->AddEntry(cJetPT,"Central Jet PT","lf");
  lg->AddEntry(fJetPT,"Forward Jet PT","lf");
  lg->Draw();

  fileName = "JetPlot-" + runNumber + ".png";
  canvas->SaveAs(fileName.c_str());

}

void MCefficiency(TString inputmc, TString outputfile, TString tfend, TString selmcgen, TString cut, TString weight, Float_t centmin, Float_t centmax)
{
  Float_t hiBinMin,hiBinMax;
  hiBinMin = centmin*2;
  hiBinMax = centmax*2;

  selmcgen = Form("%s&&hiBin>=%f&&hiBin<%f",selmcgen.Data(),hiBinMin,hiBinMax);
  cut = Form("%s&&hiBin>=%f&&hiBin<%f",cut.Data(),hiBinMin,hiBinMax);

  TFile* infMC = new TFile(inputmc.Data());
  TTree* ntMC = (TTree*)infMC->Get("ntDkpi");
  ntMC->AddFriend("ntHi");
  ntMC->AddFriend("ntSkim");
  TTree* ntGen = (TTree*)infMC->Get("ntGen");
  ntGen->AddFriend("ntHi");
  ntGen->AddFriend("ntSkim");

  TH1D* hPtMC = new TH1D("hPtMC","",nPtBins,ptBins);
  TH1D* hPtGen = new TH1D("hPtGen","",nPtBins,ptBins);

  ntMC->Project("hPtMC","Dpt",TCut(weight)*(TCut(cut.Data())&&"(Dgen==23333)"));
  divideBinWidth(hPtMC);
  ntGen->Project("hPtGen","Gpt",TCut(weight)*(TCut(selmcgen.Data())));
  divideBinWidth(hPtGen);
  TH1D* hEff = (TH1D*)hPtMC->Clone("hEff");
  hEff->Divide(hPtGen);

  TFile* fout=new TFile(Form("%s_cent_%.0f_%.0f_%s.root",outputfile.Data(),centmin,centmax,tfend.Data()),"recreate");
  fout->cd();
  hEff->Write();
  fout->Close();
}

void weightPbPbvertex(){

TFile*fMC=new TFile("/data/HeavyFlavourRun2/MC2015/Dntuple/PbPb/ntD_EvtBase_20160513_DfinderMC_PbPb_20160502_dPt1tkPt0p5_D0_prompt_Dpt2Dy1p1tkPt0p7tkEta2Decay2p9Dalpha0p14Skim_pthatweight.root");
TTree *ntDkpiMC = (TTree*)fMC->Get("ntDkpi");
TTree *ntSkimMC = (TTree*)fMC->Get("ntSkim");
TTree *ntHiMC = (TTree*)fMC->Get("ntHi");
ntDkpiMC->AddFriend(ntSkimMC);
ntDkpiMC->AddFriend(ntHiMC);

TFile*fData=new TFile("/data/jisun/PbPb2015/HF2and_ncand_skim_Dntuple_crab_PbPb_HIMinimumBias1to7_ForestAOD_D0y1p1_tkpt0p7eta1p5_goldenjson_EvtPlaneCali_03182015.root");
TTree *ntDkpiData = (TTree*)fData->Get("ntDkpi");
TTree *ntSkimData = (TTree*)fData->Get("ntSkim");
TTree *ntHiData = (TTree*)fData->Get("ntHi");
TTree *ntHltData = (TTree*)fData->Get("ntHlt");
ntDkpiData->AddFriend(ntSkimData);
ntDkpiData->AddFriend(ntHiData);
ntDkpiData->AddFriend(ntHltData);

TH1F*hpzData=new TH1F("hpzData","hpzData",200,-15,15);
TH1F*hpzMC=new TH1F("hpzMC","hpzMC",200,-15,15);

TCut weighpthat="1";
TString cut="abs(PVz)<15&&pclusterCompatibilityFilter&&pprimaryVertexFilter&&phfCoincFilter3";
TString hlt="HLT_HIL1MinimumBiasHF2AND_part1_v1||HLT_HIL1MinimumBiasHF2AND_part2_v1||HLT_HIL1MinimumBiasHF2AND_part3_v1";

TCut weighttest="1";
//TCut weighttest="1.26077+PVz*(0.00887442)+PVz*PVz*(-0.00782056)+PVz*PVz*PVz*(-3.51063e-05)+PVz*PVz*PVz*PVz*(-3.51063e-05)";

ntDkpiMC->Project("hpzMC","PVz",TCut(weighpthat)*(TCut(cut.Data()*TCut(weighttest))));
ntDkpiData->Project("hpzData","PVz",(TCut(cut.Data())*TCut(hlt.Data())));

hpzMC->Scale(1./hpzMC->Integral(hpzMC->FindBin(-15.),hpzMC->FindBin(15)));
hpzData->Scale(1./hpzData->Integral(hpzMC->FindBin(-15.),hpzMC->FindBin(15)));

TCanvas*cRatioVtx=new TCanvas("cRatioVtx","cRatioVtx",500,500);
cRatioVtx->Divide(2,1);
cRatioVtx->cd(1);
hpzMC->SetLineColor(2);
hpzMC->Draw();
hpzData->Draw("same");
cRatioVtx->cd(2);
TH1D*hRatioVertex=(TH1D*)hpzData->Clone("hRatioVertex");
hRatioVertex->Divide(hpzMC);
hRatioVertex->Draw();


TF1 *myfit = new TF1("myfit","[0]+[1]*x+x*x*[2]+x*x*x*[3]+x*x*x*x*[4]", -15, 15);  
hRatioVertex->Fit("myfit","","",-15,15);
double par0=myfit->GetParameter(0);
double par1=myfit->GetParameter(1);
double par2=myfit->GetParameter(2);
double par3=myfit->GetParameter(3);
double par4=myfit->GetParameter(4);
std::cout<<"weight="<<par0<<"+PVz*("<<par1<<")+PVz*PVz*("<<par2<<")+PVz*PVz*PVz*("<<par3<<")+PVz*PVz*PVz*PVz*("<<par4<<")"<<endl;

}

// Takes as input a given process with corresponding luminosity, as well as sample CME, cuts, directory location, and output file name - saves corresponding histogram in root file.
void processtree(TString sample, Double_t weight, TLorentzVector CME, TString filename, TString directory, TString param, TString cut, std::clock_t start, Int_t all){

    //Open File
    TFile *f = new TFile(directory);
    TTree *ttree = (TTree*)f->Get(sample);

    TString recohistname, paramhistname;
    if (all==0){
        recohistname = "_recoilmassHIST";
        paramhistname = "_parameterHIST";
    }
    if (all==1){
        recohistname = "_recoilmassHIST_all";
        paramhistname = "_parameterHIST_all";
    }

    TString histname = sample + recohistname;
    TH1D *recoilmassHIST   = new TH1D(histname,histname,CME.M()/2,0,CME.M());
    //recoilmassHIST->SetCanExtend(kAllAxes);
    ttree->Project(histname, "recoilmass", cut);

    histname = sample + paramhistname;
    TH1D *parameterHIST   = new TH1D(histname,histname,CME.M()/2,1,-1);
    //parameterHIST->SetCanExtend(kAllAxes);
    ttree->Project(histname, param, cut);

    //Normalize Hist
    recoilmassHIST->Scale(1/recoilmassHIST->GetMaximum());
    parameterHIST->Scale(1/parameterHIST->GetMaximum()) ;

    //Normalize Hist
//    recoilmassHIST->Scale(weight);
//    parameterHIST->Scale(weight) ;

    // Writes yield of each histogram for each cut to text file
    ofstream myfile;
    Double_t yield = recoilmassHIST->Integral(); 
    myfile.open ("yield.txt", ios::app);
        cout << "Yield for " << sample << " for Cut " << param << " is " << yield << endl;
        myfile << param << " " << all << " " << sample << " " << yield << "\n"; //write to file
    myfile.close();

    if (recoilmassHIST->Integral() == 0){
        recoilmassHIST->Fill(0);
    }
    if (parameterHIST->Integral() == 0){
        parameterHIST->Fill(0);
    }

    TFile g(filename, "update");
    recoilmassHIST->Write();
    parameterHIST->Write();
    g.Close();
}

void plots_for_config(string toyfile, string smsfile, string massGlu, string massLSP) {

  gROOT->SetStyle("Plain");
  bool do_fat_line=false; // if you want to have HistLineWidth=1 and FuncWidth=1 as it was before instead of 2
  setTDRStyle(do_fat_line);
  gStyle->SetTextFont(42);
  
  string massCut_ = "MassGlu==" + massGlu + "&&MassLSP==" + massLSP;
  string massCut2_ = "massGlu==" + massGlu + "&&massLSP==" + massLSP;
  TCut massCut(massCut_.c_str());
  TCut massCut2(massCut2_.c_str());

  string comp = "comparison_" + massGlu + "_" + massLSP;

  TFile *_file = TFile::Open(smsfile.c_str());
  TTree *point = (TTree*)_file->Get("events");
  TH1F *jzb = new TH1F("jzb", "", 50, -600, 800);
  point->Project("jzb", "jzb[1]", "id1==id2&&pfJetGoodNum>2&&abs(mll-91.2)<20" + massCut, "gof");

  TFile *_filepa = TFile::Open(toyfile.c_str());
  TTree *pointpa = (TTree*)_filepa->Get("events");
  TH1F *jzbpa = new TH1F("jzbpa", "", 50, -600, 800);
  pointpa->Project("jzbpa", "jzb", "njets>2 && nleptons==2" + massCut2);


  jzbpa->SetLineColor(kRed);
  jzbpa->GetXaxis()->SetTitle("JZB [GeV]");
  jzb->GetXaxis()->SetTitle("JZB [GeV]");


  TLegend *l1 = new TLegend(0.7, 0.75, 0.9, 0.85);
  l1->SetFillColor(kWhite);
  l1->SetTextFont(42);
  l1->SetLineWidth(0);
  l1->SetBorderSize(0);
  l1->AddEntry(jzb, "SMS", "L");
  l1->AddEntry(jzbpa, "SMS TOY", "L");
  TCanvas *m1 = new TCanvas("m1");
  m1->GetPad(0)->SetLogy(1);
  m1->cd();
  jzb->DrawNormalized();
  jzbpa->DrawNormalized("SAME");
  l1->Draw(); 
  CompleteSave(m1,"Plots/" + comp);
  delete m1;

  return;

}

void ROCTnP(){
  gROOT->LoadMacro("tdrstyle.C");
  setTDRStyle();

  TFile * fbkg = new TFile("/home/tjkim/ntuple/iso/data/QCD_tight/vallot.root");
  TTree * tbkg = (TTree *) fbkg->Get("MuonAna/tree");

  //from normalized histos
  TFile * fhisto = new TFile("MuonAna_backup.root");
  TH1F * h_TTbar = (TH1F *) fhisto->Get("Step_1/hMCSig_TTbar_Step_1_relIso");
  TH1F * h_Wm = (TH1F *) fhisto->Get("Step_1/hMC_Wm_Step_1_relIso");
  TH1F * h_Zmm = (TH1F *) fhisto->Get("Step_1/hMC_Zmumu_Step_1_relIso");
  TH1F * h_QCD = (TH1F *) fhisto->Get("Step_1/hMC_QCD_Step_1_relIso");
  TH1F * h_Data = (TH1F *) fhisto->Get("Step_1/hData_Step_1_relIso");

  h_QCD->Add(h_TTbar);  
  h_QCD->Add(h_Wm);  
  h_QCD->Add(h_Zmm);  

  TCut cutbkg = "pt > 20 && pt < 50"; 
  
  TH1* h_bkg = new TH1F("h_bkg","h_bkg",50,0,0.5);
  tbkg->Project("h_bkg","(chIsoOpt04+phIsoOpt04+nhIsoOpt04)/pt", cutbkg);

  //Data vs PU
  plot(h_Data, h_QCD, "Data (T&P)", "MC (T&P)" , "QCD_MC_pt20_50_PU" , "p_{T}=20-50 GeV" );
  
}

void leadFrag(char *infname = "dj_HCPR-J50U-hiGoodMergedTracks_OfficialSelv2_Final0.root",
    TCut myCut = "cent<30",
    TString module="djcalo_tower")
{
   TFile *inf = new TFile(infname);
   TTree *t = (TTree*)inf->Get(module+"/djTree");

   const int nBin = 5;
   double bins[nBin+1] = {0.5,1.5,4,8,20,1000};  
   double colors[nBin] = {kGray,kBlue-3,38,kOrange-8,kRed-6};

   // Selection cut
   TCut evtCut = "nljet>120&&abs(nljeta)<2&&aljet>50&&abs(aljeta)<2&&jdphi>2./3*TMath::Pi()";
   t->SetAlias("Aj","(nljet-aljet)/(nljet+aljet)");

   // histograms
   TH2D * hFrag_LJEt = new TH2D("hFrag_LJEt",";p_{T}^{J1};p_{T}^{Tower}/p_{T}^{J1}",50,120,300,50,0,1);
   hFrag_LJEt->GetXaxis()->CenterTitle();
   hFrag_LJEt->GetYaxis()->CenterTitle();
   TH2D * hLeadFrag_LJEt = new TH2D("hLeadFrag_LJEt",";p_{T}^{J1};p_{T}^{Highest Tower}/p_{T}^{J1}",50,120,300,50,0,1);
   hLeadFrag_LJEt->GetXaxis()->CenterTitle();
   hLeadFrag_LJEt->GetYaxis()->CenterTitle();

   // project tree
   t->Project("hLeadFrag_LJEt","lppt[0]/nljet:nljet",evtCut);
   t->Project("hFrag_LJEt","ppt/nljet:nljet",evtCut);


   // plot
   TCanvas *c2 = new TCanvas("c2","",600,600);
   c2->SetRightMargin(0.15);
   c2->SetLogz();
   hFrag_LJEt->Draw("colz");

   TCanvas *c3 = new TCanvas("c3","",600,600);
   c3->SetRightMargin(0.15);
   c3->SetLogz();
   hLeadFrag_LJEt->Draw("colz");

   TCanvas *c4 = new TCanvas("c4","",600,600);
   hLeadFrag_LJEt->ProfileX()->Draw("E");
   TH1D * hLeadFrag_LJEt_pfx = (TH1D*)gDirectory->Get("hLeadFrag_LJEt_pfx");
   hLeadFrag_LJEt_pfx->SetAxisRange(0,0.6,"Y");
   hFrag_LJEt->ProfileX()->Draw("same hist");
   TH1D * hFrag_LJEt_pfx = (TH1D*)gDirectory->Get("hFrag_LJEt_pfx");
}

void get_count_and_uncertainty(TTree& tree,
                               const string& cut,
                               double& count,
                               double& uncertainty) {
    const string hist_name("temp");
    TH1D temp(hist_name.c_str(), "", 1, -1.0, 1.0);
    tree.Project(hist_name.c_str(), "0.0", cut.c_str());
    count=temp.IntegralAndError(0,2,uncertainty);
}

void addbcsv(){

  gROOT->ProcessLine(".L $CMSSW_BASE/src/KoPFA/CommonTools/macros/tdrstyle.C");
  defaultStyle();

  TFile * f = new TFile("/afs/cern.ch/work/y/youngjo/public/For8Tev/v20130612_genjet/vallot_TTbarFullLepMGDecays.root");
  TTree * t = (TTree*) f->Get("MuEl/tree");
  
  TH1F * h_bfromtop =  new TH1F("h_bfromtop","h_bfromtop",20,0,1);
  TH1F * h_bnotfromtop =  new TH1F("h_bnotfromtop","h_bnotfromtop",20,0,1);

  t->Project("h_bfromtop","jets_bDiscriminatorCSV","jets_fromtop == 1 && abs(jets_flavor) == 5");
  t->Project("h_bnotfromtop","jets_bDiscriminatorCSV","jets_fromtop == 0 && abs(jets_flavor) == 5");

  h_bfromtop->Scale(1.0/h_bfromtop->Integral());
  h_bnotfromtop->Scale(1.0/h_bnotfromtop->Integral());

  TCanvas * c = new TCanvas("c","c",1);

  h_bfromtop->Draw();
  h_bfromtop->SetStats(0);
  h_bfromtop->SetTitle("");
  h_bfromtop->SetLineColor(2);
  h_bfromtop->SetLineWidth(2);
  h_bnotfromtop->Draw("same");
  h_bnotfromtop->SetLineColor(4);
  h_bnotfromtop->SetLineWidth(2);

  TLegend *l = new TLegend(0.48,0.65,0.58,0.80);
  l->AddEntry(h_bfromtop,"b from top","L");
  l->AddEntry(h_bnotfromtop,"b not from top","L");
  l->SetTextSize(0.04);
  l->SetFillColor(0);
  l->SetLineColor(0);

  l->Draw();

  TLatex *label= new TLatex;
  label->SetNDC();
  label->SetTextSize(0.05);
  label->DrawLatex(0.20,0.82,"CMS Simulation, at #sqrt{s} = 8 TeV");

  c->Print("csvshape_addb.eps");
}

void GetCountAndUncertainty(TTree &tree,
                            const Cut &cut,
                            double &count,
                            double &uncertainty){
  const string hist_name{"temp"};
  TH1D temp{hist_name.c_str(), "", 1, -1.0, 1.0};
  temp.Sumw2();
  tree.Project(hist_name.c_str(), "0.", static_cast<const char *>(cut));
  count=temp.IntegralAndError(0,2,uncertainty);
}

void photonSpectra()
{
  TH1::SetDefaultSumw2();

  const Double_t PTBINS[] = {40, 50, 60, 80, 120, 300};
  const Int_t nPTBINS = 5;

  const TString DATA_FILE = "gammaJets_pA_Data.root";
  const TString MC_FILE = "gammaJets_pA_MC_allQCDPhoton.root";
  TFile *dataFile = TFile::Open(DATA_FILE);
  TTree *dataTree = (TTree*)dataFile->Get("photonTree");
  TFile *mcFile = TFile::Open(MC_FILE);
  TTree *mcTree = (TTree*)mcFile->Get("photonTree");

  TH1D *photonSpectra = new TH1D("photonSpectra","photonSpectra", nPTBINS, PTBINS);
  TH1D *correctedPhotonSpectra = (TH1D*)photonSpectra->Clone("correctedPhotonSpectra");

  const TCut photonEtaCut = "abs(photonTree.eta) < 1.44";
  const TCut isoCut = "ecalRecHitSumEtConeDR04 < 4.2  &&  hcalTowerSumEtConeDR04 < 2.2  &&  trkSumPtHollowConeDR04 < 2 && hadronicOverEm<0.1";
  const TCut genCut = "genCalIsoDR04<5 && abs(genMomId)<=22";
  const TCut sbCut = "(cc4+cr4+ct4PtCut20>10) && (cc4+cr4+ct4PtCut20<20) && hadronicOverEm<0.1";
  const TCut candidateCut = "sigmaIetaIeta<0.01";
  const TCut decayCut = "(sigmaIetaIeta>0.011) && (sigmaIetaIeta<0.017)";

  TCut photonPtCut = "photonTree.corrPt>40 && photonTree.corrPt<300";

  dataTree->Project(photonSpectra->GetName(),"corrPt",isoCut && photonEtaCut && photonPtCut && candidateCut);

  for(int i = 0; i<nPTBINS; ++i)
  {
    photonPtCut = Form("photonTree.corrPt>%f && photonTree.corrPt<%f", PTBINS[i], PTBINS[i+1]  );

    TCut dataCandidateCut = isoCut && photonEtaCut && photonPtCut;
    TCut sidebandCut =  sbCut && photonEtaCut && photonPtCut;
    TCut mcSignalCut = dataCandidateCut && genCut;

    fitResult fitr = getPurity(dataTree, mcTree,
  			       dataCandidateCut, sidebandCut,
  			       mcSignalCut);
    Double_t purity = fitr.purity;
    correctedPhotonSpectra->SetBinContent(i+1,photonSpectra->GetBinContent(i+1)*purity);
  }

  // TCanvas *c1 = new TCanvas();
  // photonSpectra->Draw();
  // correctedPhotonSpectra->Draw("same");

  TFile *outFile = new TFile("photonSpectra.root", "RECREATE");
  outFile->cd();
  photonSpectra->Write();
  correctedPhotonSpectra->Write();
  outFile->Close();
}

void cdcHits(void)
{
	gROOT->Reset();
	gStyle->SetErrorX(0.0);

	TCanvas *c1 = new TCanvas("c1");
	c1->SetTicky();
	c1->SetTickx();
	c1->SetGridy();
	c1->Draw();
	
	TFile *f = new TFile("cdc_hits.root");
	TTree *t = (TTree*)gROOT->FindObject("cdcHits");
	double Nevents_thrown = (double)t->GetEntries();
	double timewindow = 1.0E-6;
	double total_integrated_time = timewindow*Nevents_thrown;
	double tagged_rate = 1.0E7;
	cout<<"Number of events thrown: "<<(int)Nevents_thrown<<endl;
	cout<<"Assuming "<<timewindow/1.0E-6<<" microsecond time window"<<endl;
	cout<<"Total beam time represented by simulation:"<<total_integrated_time<<" seconds"<<endl;
	cout<<"Assuming rates correspond to "<<tagged_rate<<" tagged photons/sec beam"<<endl;
	
	// Number of straws in each "ring"
	int n_straws[]={43,50,57,64,71,78,85,99,106,113,120,127,134,141,148,155,166,173,182,187,194,201,208,215,222};
	
	TH1D *nstraws = new TH1D("nstraws","Straws per ring", 25, 0.5, 25.5);
	for(int i=0;i<25; i++){
		nstraws->SetBinContent(i+1, n_straws[i]);
		nstraws->SetBinError(i+1, 0.0);
	}

	TH1D *hits_per_straw = new TH1D("hits_per_straw","CDC Hits per layer", 25, 0.5, 25.5);
	t->Project("hits_per_straw","ring");
	hits_per_straw->Divide(nstraws);
	
	TH1D *rate_per_straw = hits_per_straw->Clone("rate_per_straw");
	char str[256];
	sprintf(str, "Rate per straw at %g tags/sec (kHz)", tagged_rate);
	rate_per_straw->Scale(1.0E-3/total_integrated_time);
	rate_per_straw->SetXTitle("CDC layer");
	rate_per_straw->SetYTitle(str);
	rate_per_straw->SetStats(0);
	rate_per_straw->SetLineColor(kRed);
	rate_per_straw->SetLineWidth(2.0);
	rate_per_straw->SetMarkerColor(kRed);
	rate_per_straw->SetMarkerStyle(20);
	rate_per_straw->SetMarkerSize(1.5);
	rate_per_straw->Draw("P");
	
	c1->SaveAs("cdc_hit_rate.gif");
	c1->SaveAs("cdc_hit_rate.pdf");
}

void saveHist3pMinpt()
{
    real isMC=MC;
    TFile* infMCData = new TFile(infnameData3p[isMC].Data());
    TFile* infMCMC = new TFile(infnameMC3p[isMC].Data());
    TTree* ntMCData = (TTree*)infMCData->Get("ntDD0kpipi");
    TTree* ntMCMC = (TTree*)infMCMC->Get("ntDD0kpipi");

    ntMCData->AddFriend("ntHlt");
    if(isMC!=Data_MB||isMC!=Data) ntMCData->AddFriend("ntHi");
    ntMCMC->AddFriend("ntHlt");
    ntMCMC->AddFriend("ntHi");

    cout<<"  -- Filling histograms - MC"<<endl;
    cout<<"     "<<infnameData3p[isMC]<<endl;
    cout<<"     "<<infnameMC3p[isMC]<<endl;

    for(int i=0; i<nBins; i++)
    {
        cout<<setiosflags(ios::left)<<"   - Processing ptbin "<<setw(3)<<ptBins[i]<<endl;
        TH1D* h = new TH1D("h","",BINNUM,BINMIN,BINMAX);
        TH1D* hMCSignal = new TH1D("hMCSignal","",BINNUM,BINMIN,BINMAX);
        TH1D* hMCSwapped = new TH1D("hMCSwapped","",BINNUM,BINMIN,BINMAX);
        ntMCData->Project("h","Dmass-DtktkResmass",Form("%s*(%s&&%s&&Dpt>%f)",weightdata[isMC].Data(),seldata3p[isMC].Data(),triggerselectiondata[isMC].Data(),ptBins[i]));
        ntMCMC->Project("hMCSignal","Dmass-DtktkResmass",Form("%s*(%s&&%s&&Dpt>%f)",weightmc[isMC].Data(),selmc3p[isMC].Data(),triggerselectionmc[isMC].Data(),ptBins[i]));
        ntMCMC->Project("hMCSwapped","Dmass-DtktkResmass",Form("%s*(%s&&%s&&Dpt>%f)",weightmc[isMC].Data(),selswp3p[isMC].Data(),triggerselectionmc[isMC].Data(),ptBins[i]));
        TFile* outf = new TFile(Form("inputfiles/fmass_pp_3p_%s_%.0f.root",texData[isMC].Data(),ptBins[i]),"recreate");
        outf->cd();
        h->Write();
        hMCSignal->Write();
        hMCSwapped->Write();
        outf->Close();
        delete hMCSwapped;
        delete hMCSignal;
        delete h;
        delete outf;
    }
    cout<<endl;
}

void timeWindow( TFile *infile, double w0, double w1, double w2, double w3)
{
	// w0 -> w1 = time window around peak of time-expt
	// w2 -> w0 = window to left of w01
	// w1 -> w3 = window to right of w01

	// grab trees from infile
	TTree *dirc = infile->Get("dirc");
	TTree *reco = infile->Get("reco");

	// get lens, bar angle, and beam momentum
	int lens(0);
	double alpha(0),beam(0);
	dirc->SetBranchAddress("lens",&lens);
	dirc->SetBranchAddress("theta",&alpha);
	dirc->GetEntry(0);

	// define histograms of diff, outside diff, and counters
	TH1D *diff = new TH1D("diff","diff",1000,-50,50);
	TH1D *dout = new TH1D("dout","dout",1000,-50,50);
	TH1D *tof1 = new TH1D("time-tof1","time-tof1",1000,0,1000);
	TH1D *tof2 = new TH1D("time-tof2","time-tof2",1000,0,1000);
	TH1D *trig = new TH1D("time-trig","time-trig",1000,0,1000);

	// project to histograms
	reco->Project("diff","diff");
	reco->Project("dout","diff");
	reco->Project("time-tof1","time-tof1");
	reco->Project("time-tof2","time-tof2");
	reco->Project("time-trig","time-trig");

	// get diff value at peak and define initial cut
	double diffpeak = diff->GetMaximumBin();
	diffpeak = (diffpeak-500)/10;
	char *diffcut = Form("%f<diff-%f && diff-%f<%f",w0,diffpeak,diffpeak,w1);
	cout << "Initial cut\t" << diffcut << endl;
	
	
}