Plotting/Modules/FSRStudy.C

#include <iostream>


samplecollection PhotonSamples("PhotonSamples");

namespace GenPhotons {
  void CompareZmassToFinalStateLeptons();
  void GenLevelStudies();
  
  void NumberOfPhotonsPerEvent();
  void CharacterizingFSRPhotons();
  void MakeCharacterizationPlot(string variable, int nbins, float min, float max, string label, string saveas);
  void dR2DPlots();
  void CompareGenJZBDistributions();
  void EfficiencyPurity();
    
  TCut FSRcut;
  TCut NoFSRcut;
  TCut BaseCut;
}

namespace RecoPhotons {
  void RecoLevelStudies();
  
  
}

void RecoPhotons::RecoLevelStudies() {
  cout << "RecoLevelStudies have not been implemented yet" << endl;
}


void GenPhotons::CompareZmassToFinalStateLeptons() {
  
  TCanvas *can = new TCanvas("can","can");
  can->SetLogy(1);
  
  TH1F *FullZ = PhotonSamples.Draw("FullZ","genFSRmllg",36,20,200,"m_{Z}^{gen} [GeV]",  "events", "genFSRmllg>20&&genFSRmllg<200",mc,PlottingSetup::luminosity);
  TH1F *GenLL = PhotonSamples.Draw("GenLL","genFSRmll", 36,20,200,"m_{ll}^{gen} [GeV]", "events","genFSRmllg>20&&genFSRmllg<200",mc,PlottingSetup::luminosity);
  
  FullZ->SetLineColor(TColor::GetColor("#01DF01"));
  GenLL->SetLineColor(TColor::GetColor("#FF4000"));
  
  TLegend *leg = new TLegend(0.7,0.7,0.89,0.89);
  leg->SetFillColor(kWhite);
  leg->AddEntry(FullZ,"gen Z mass","l");
  leg->AddEntry(GenLL,"gen m_{ll}","l");
  
  FullZ->Draw("histo");
  GenLL->Draw("same,histo");
  leg->SetHeader("DY MC:");
  leg->Draw("same,histo");
  DrawMCPrelim();
  
  CompleteSave(can,"PhotonStudies/CompareZmassToFinalStateLeptons");
  delete can;
}  


void GenPhotons::NumberOfPhotonsPerEvent() {
    TCanvas *can = new TCanvas("can","can");
    can->SetLogy(1);
    
    TH1F *FSR = PhotonSamples.Draw("FSR","genFSRNparticles"  ,11,-0.5,10.5,"N(#gamma)", "events", BaseCut,mc,PlottingSetup::luminosity);
    TH1F *all = PhotonSamples.Draw("all","genPhotonsNPhotons",11,-0.5,10.5,"N(#gamma)", "events", BaseCut,mc,PlottingSetup::luminosity);
    
    FSR->SetLineColor(TColor::GetColor("#01DF01"));
    all->SetLineColor(TColor::GetColor("#FF4000"));
    
    TLegend *leg = new TLegend(0.7,0.7,0.89,0.89);
    leg->SetFillColor(kWhite);
    leg->AddEntry(FSR,"FSR photons","l");
    leg->AddEntry(all,"all photons","l");
    
    FSR->Draw("histo");
    all->Draw("same,histo");
    leg->SetHeader("DY MC:");
    leg->Draw("same,histo");
    DrawMCPrelim();
    
    CompleteSave(can,"PhotonStudies/NPhotons");
    
    
    TH1F *PhotonPt[10];
    for(int i=0;i<10;i++) {
        PhotonPt[i] = PhotonSamples.Draw("FSR","genPhotonsPt",25,0,100,"N(#gamma)",  "events", FSRcut && TCut(("genFSRNparticles=="+any2string(i)).c_str()),mc,PlottingSetup::luminosity);
    }
    
    int DrawCounter=0;
    
    TLegend *leg2 = make_legend();
    for(int i=0;i<10;i++) {
        cout << PhotonPt[i]->Integral() << endl;
        if(PhotonPt[i]->Integral()>50) {
            PhotonPt[i]->SetLineColor(DrawCounter+1);
            leg2->AddEntry(PhotonPt[i],(any2string(i)+" FSR Photons").c_str(),"l");
            if(DrawCounter) PhotonPt[i]->DrawNormalized("same,hist");
            else PhotonPt[i]->DrawNormalized("hist");
            DrawCounter++;
        }
    }
    
    leg2->Draw();
    
    leg2->SetHeader("DY MC (N_{jets} #geq 0):");
    DrawMCPrelim();
    
    CompleteSave(can,"PhotonStudies/PhotonPtvsNPhotons");
    
    delete can;
}

void GenPhotons::MakeCharacterizationPlot(string variable, int nbins, float min, float max, string label, string saveas) {
  TCanvas *can = new TCanvas("can","can");
  can->SetLogy(1);
  
  TH1F *dR_FSR  = PhotonSamples.Draw("dR_FSR", variable, nbins, min, max, label, "events",BaseCut && FSRcut,  mc,PlottingSetup::luminosity);
  TH1F *dR_nFSR = PhotonSamples.Draw("dR_nFSR", variable, nbins, min, max, label, "events",BaseCut && NoFSRcut, mc,PlottingSetup::luminosity);
  
  dR_FSR->SetLineColor(TColor::GetColor("#01DF01"));
  dR_nFSR->SetLineColor(TColor::GetColor("#FF4000"));
  
  TLegend *leg = new TLegend(0.7,0.7,0.89,0.89);
  leg->SetFillColor(kWhite);
  leg->AddEntry(dR_FSR,"FSR #gamma's","l");
  leg->AddEntry(dR_nFSR,"non-FSR #gamma's","l");
  
  if(dR_nFSR->GetMinimum()<dR_FSR->GetMinimum()) dR_FSR->SetMinimum(dR_nFSR->GetMinimum());
  if(dR_nFSR->GetMaximum()>dR_FSR->GetMaximum()) dR_FSR->SetMaximum(dR_nFSR->GetMaximum());
  
  
  dR_FSR->Draw("histo");
  dR_nFSR->Draw("same,histo");
  leg->SetHeader("DY MC:");
  leg->Draw("same,histo");
  DrawMCPrelim();
  
  CompleteSave(can,"PhotonStudies/CharacterizingFSRPhotons/"+saveas);
  delete dR_FSR;
  delete dR_nFSR;
  delete can;
}


void GenPhotons::CharacterizingFSRPhotons() {
  MakeCharacterizationPlot("min(genPhotonsdR1,genPhotonsdR2)",100,0,5,"min(#Delta R (#gamma,l_{1}), #Delta R (#gamma,l_{2}))", "MinDR");
  MakeCharacterizationPlot("genPhotonsPt",100,0,100,"p_{T}^{#gamma}", "PhotonPt");
  MakeCharacterizationPlot("genPhotonsEta",100,-5,5,"{#eta}^{#gamma}", "PhotonEta");
  MakeCharacterizationPlot("genPhotonsPhi",100,0,4,"{#phi}^{#gamma}", "PhotonPhi");
  MakeCharacterizationPlot("genPhotonsM",100,0,2,"m^{#gamma}", "PhotonM");
}

void GenPhotons::CompareGenJZBDistributions() {
    TCanvas *can = new TCanvas("can","can");
    can->SetLogy(1);
    
    TH1F *dR_FSR  = PhotonSamples.Draw("dR_FSR", "genFSRjzbG", 30, -30, 30, "JZB^{gen}", "events",BaseCut && TCut("abs(genFSRjzbG)>0.1"),   mc,PlottingSetup::luminosity);
    TH1F *dR_nFSR = PhotonSamples.Draw("dR_nFSR", "genFSRjzb", 30, -30, 30, "JZB^{gen}", "events",BaseCut && TCut("abs(genFSRjzb)>0.1") , mc,PlottingSetup::luminosity);
    
    dR_FSR->SetLineColor(TColor::GetColor("#01DF01"));
    dR_nFSR->SetLineColor(TColor::GetColor("#FF4000"));
    
    TLegend *leg = new TLegend(0.7,0.7,0.89,0.89);
    leg->SetFillColor(kWhite);
    leg->AddEntry(dR_FSR,"With FSR recovery","l");
    leg->AddEntry(dR_nFSR,"Without FSR recovery","l");
    
    if(dR_nFSR->GetMinimum()<dR_FSR->GetMinimum()) dR_FSR->SetMinimum(dR_nFSR->GetMinimum());
    if(dR_nFSR->GetMaximum()>dR_FSR->GetMaximum()) dR_FSR->SetMaximum(dR_nFSR->GetMaximum());
    
    
    dR_FSR->Draw("histo");
    dR_nFSR->Draw("same,histo");
    leg->SetHeader("DY MC:");
    leg->Draw("same,histo");
    DrawMCPrelim();
    
    CompleteSave(can,"PhotonStudies/GenJZBDistributions");
    delete dR_FSR;
    delete dR_nFSR;
    delete can;
}

void GenPhotons::dR2DPlots() {
    TCanvas *can = new TCanvas("can","can",1000,500);
    can->Divide(2,1);
    can->cd(1);
    
    vector<float> binX, binY;
    for(int i=0;i<=31;i++) {
        binX.push_back(i*0.1);
        binY.push_back(i*0.1);
    }
    
    TH2F *dR_FSR   = PhotonSamples.Draw("dR_FSR",  "genPhotonsdR1:genPhotonsdR2", binX,binY, "#Delta R (#gamma^{FSR},l_{1})","#Delta R (#gamma^{FSR},l_{2})",BaseCut && FSRcut   ,mc,PlottingSetup::luminosity, PhotonSamples.FindSample("o"),false);
    TH2F *dR_nFSR  = PhotonSamples.Draw("dR_nFSR", "genPhotonsdR1:genPhotonsdR2", binX,binY, "#Delta R (#gamma^{non-FSR},l_{1})","#Delta R (#gamma^{non-FSR},l_{2})",BaseCut && NoFSRcut ,mc,PlottingSetup::luminosity, PhotonSamples.FindSample("o"),false);
    
    can->cd(1);
    dR_FSR->Draw("COLZ");
    DrawMCPrelim();
    
    can->cd(2);
    dR_nFSR->Draw("COLZ");
    DrawMCPrelim();
    
    CompleteSave(can,"PhotonStudies/dR_in_2D");
    delete dR_FSR;
    delete dR_nFSR;
    delete can;
}


void GenPhotons::EfficiencyPurity() {
    TCanvas *can = new TCanvas("can","can");
    cout << "Computing efficiency and purity !" << endl;
    vector<float> bindR, binPt;
    for(int i=0;i<=100;i++) bindR.push_back(i*0.05);
    binPt.push_back(0);binPt.push_back(10);binPt.push_back(50);binPt.push_back(100);binPt.push_back(1000);binPt.push_back(1010);
//    binPt.push_back(0);binPt.push_back(1000);binPt.push_back(1001);
    
    
    TGraph *purity[binPt.size()];
    TGraph *efficiency[binPt.size()];
    TGraph *purity_n_efficiency[binPt.size()];
    TH1F *FSR[binPt.size()];
    TH1F *Any[binPt.size()];
    for(int i=1;i<binPt.size();i++) {
       //do this for each pt range
       float ptlow=binPt[i-1];
       float pthi =binPt[i];
       stringstream ptcut;
       ptcut << "genPhotonsPt>" << ptlow << " && genPhotonsPt<" << pthi << " && min(genPhotonsdR1,genPhotonsdR2)<5";
       
       FSR[i-1]  = PhotonSamples.Draw("FSR", "min(genPhotonsdR1,genPhotonsdR2)", bindR,"min DR","events",TCut(ptcut.str().c_str()) && BaseCut && FSRcut ,mc,PlottingSetup::luminosity);
       Any[i-1]  = PhotonSamples.Draw("Any", "min(genPhotonsdR1,genPhotonsdR2)", bindR,"min DR","events",TCut(ptcut.str().c_str()) && BaseCut  ,mc,PlottingSetup::luminosity);
       
       float NPhotons=0;
       float NFSRPhotons=0;
       float totalNFSRphotons=0;
       
       for(int k=1;k<FSR[i-1]->GetNbinsX();k++) totalNFSRphotons+=FSR[i-1]->GetBinContent(k);
       cout << i << " : " << totalNFSRphotons << " FSR photons in pt range " << ptlow << " , " << pthi << endl;
       purity[i-1]=new TGraph(FSR[i-1]->GetNbinsX());
       efficiency[i-1]=new TGraph(FSR[i-1]->GetNbinsX());
       purity_n_efficiency[i-1]=new TGraph(FSR[i-1]->GetNbinsX());
       for(int j=1;j<=FSR[i-1]->GetNbinsX();j++) {
         NPhotons+=(int)Any[i-1]->GetBinContent(j);
         NFSRPhotons+=FSR[i-1]->GetBinContent(j);
//            cout << "       ::::::::" << j << " : " << NPhotons << " : " << NFSRPhotons << " : " << totalNFSRphotons << endl;
            if(NPhotons>0) purity[i-1]->SetPoint(j-1,FSR[i-1]->GetXaxis()->GetBinCenter(j)+FSR[i-1]->GetBinWidth(j),(float)NFSRPhotons/NPhotons);
            else purity[i-1]->SetPoint(j-1,FSR[i-1]->GetXaxis()->GetBinCenter(j)+FSR[i-1]->GetBinWidth(j),-1);
            if(totalNFSRphotons>0) efficiency[i-1]->SetPoint(j,FSR[i-1]->GetXaxis()->GetBinCenter(j)+FSR[i-1]->GetBinWidth(j),(float)NFSRPhotons/totalNFSRphotons);
            else efficiency[i-1]->SetPoint(j-1,FSR[i-1]->GetXaxis()->GetBinCenter(j)+FSR[i-1]->GetBinWidth(j),-1);
            if(totalNFSRphotons>0&&NPhotons>0) purity_n_efficiency[i-1]->SetPoint(j-1,(float)NFSRPhotons/totalNFSRphotons,(float)NFSRPhotons/NPhotons);
        }
        cout << "Done for pt range" << endl;
    }
    cout << "Done with graphs" << endl;
//    cout << "Artificially aborting ... " << endl;return;
    for(int i=0;i<binPt.size()-1;i++) {
      cout << "About to draw purity for i=" << i << endl;
      cout << "( this corresponds to the Pt range from " << binPt[i] << " to " << binPt[i+1] << ")" << endl;
        if(i==0) purity[i]->Draw("AC");
        else purity[i]->Draw("C");
    }
    CompleteSave(can,"PhotonStudies/Purity");
    
    for(int i=0;i<binPt.size()-1;i++) {
      cout << "About to draw efficiency for i=" << i << endl;
        if(i==0) efficiency[i]->Draw("AC");
        else efficiency[i]->Draw("C");
    }
    CompleteSave(can,"PhotonStudies/Efficiency");
    
    for(int i=0;i<binPt.size()-1;i++) {
        if(i==0) purity_n_efficiency[i]->Draw("AC");
        else purity_n_efficiency[i]->Draw("C");
    }
    CompleteSave(can,"PhotonStudies/Purity_n_efficiency");
}


void GenPhotons::GenLevelStudies() {
    TCut essential_bkp = essentialcut;
    essentialcut=TCut("");
    
//  CompareZmassToFinalStateLeptons();
//  NumberOfPhotonsPerEvent();
//  CharacterizingFSRPhotons();
//  CompareGenJZBDistributions();
//  dR2DPlots();
    
  EfficiencyPurity();
  
  essentialcut=essential_bkp;
}

void InitializePhotonSamples() {
  float ZJetsCrossSection         = 3532.8;
  float LowZJetsCrossSection      = 11050*0.069*1.15; // LO xs * filter eff * k-factor (same as ZJets)

  PhotonSamples.AddSample("/shome/buchmann/ntuples/MC8tev/PhotonStudies/Skim3genJets_DYJetsToLL_M-50_TuneZ2Star_8TeV-madgraph-tarball-Summer12_DR53X-PU_S10_START53_V7A-v1.root","Z+Jets",-30346766,ZJetsCrossSection,false,false,1,kYellow);
  PhotonSamples.AddSample("/shome/buchmann/ntuples/MC8tev/PhotonStudies/Skim3genJets_DYJetsToLL_M-10To50filter_8TeV-madgraph-Summer12_DR53X-PU_S10_START53_V7A-v1.root","Z+Jets",-6955271,LowZJetsCrossSection,false,false,1,kYellow);
  
  write_info(__FUNCTION__,"Photon samples have been initialized!");
  PhotonSamples.ListSamples();
  
  GenPhotons::FSRcut   = TCut("genPhotonsIsFSR");
  GenPhotons::NoFSRcut = TCut("!genPhotonsIsFSR");
  GenPhotons::BaseCut = TCut("genNjets>=3 && genPhotonsPt > 5 && NgenLeps>0");
}

void FSRstudy() {
  InitializePhotonSamples();
  
  GenPhotons::GenLevelStudies();
  RecoPhotons::RecoLevelStudies();
}
Revision:	1.3
Committed:	Mon Oct 15 22:12:14 2012 UTC (12 years, 6 months ago) by buchmann
Content type:	text/plain
Branch:	MAIN
Changes since 1.2:	+45 -48 lines
Log Message:	Updated efficiency/purity plots: a bit slower now but much easier to read and adapt
#	User	Rev	Content
1	buchmann	1.1	#include <iostream>
2
3
4			samplecollection PhotonSamples("PhotonSamples");
5
6			namespace GenPhotons {
7			void CompareZmassToFinalStateLeptons();
8			void GenLevelStudies();
9
10	buchmann	1.2	void NumberOfPhotonsPerEvent();
11			void CharacterizingFSRPhotons();
12			void MakeCharacterizationPlot(string variable, int nbins, float min, float max, string label, string saveas);
13			void dR2DPlots();
14			void CompareGenJZBDistributions();
15			void EfficiencyPurity();
16
17			TCut FSRcut;
18			TCut NoFSRcut;
19			TCut BaseCut;
20	buchmann	1.1	}
21
22			namespace RecoPhotons {
23			void RecoLevelStudies();
24
25
26
27			}
28
29			void RecoPhotons::RecoLevelStudies() {
30			cout << "RecoLevelStudies have not been implemented yet" << endl;
31			}
32
33
34			void GenPhotons::CompareZmassToFinalStateLeptons() {
35
36			TCanvas *can = new TCanvas("can","can");
37			can->SetLogy(1);
38
39	buchmann	1.2	TH1F *FullZ = PhotonSamples.Draw("FullZ","genFSRmllg",36,20,200,"m_{Z}^{gen} [GeV]", "events", "genFSRmllg>20&&genFSRmllg<200",mc,PlottingSetup::luminosity);
40			TH1F *GenLL = PhotonSamples.Draw("GenLL","genFSRmll", 36,20,200,"m_{ll}^{gen} [GeV]", "events","genFSRmllg>20&&genFSRmllg<200",mc,PlottingSetup::luminosity);
41	buchmann	1.1
42			FullZ->SetLineColor(TColor::GetColor("#01DF01"));
43			GenLL->SetLineColor(TColor::GetColor("#FF4000"));
44
45			TLegend *leg = new TLegend(0.7,0.7,0.89,0.89);
46	buchmann	1.2	leg->SetFillColor(kWhite);
47	buchmann	1.1	leg->AddEntry(FullZ,"gen Z mass","l");
48	buchmann	1.2	leg->AddEntry(GenLL,"gen m_{ll}","l");
49	buchmann	1.1
50	buchmann	1.2	FullZ->Draw("histo");
51			GenLL->Draw("same,histo");
52			leg->SetHeader("DY MC:");
53			leg->Draw("same,histo");
54			DrawMCPrelim();
55	buchmann	1.1
56			CompleteSave(can,"PhotonStudies/CompareZmassToFinalStateLeptons");
57	buchmann	1.2	delete can;
58			}
59
60
61			void GenPhotons::NumberOfPhotonsPerEvent() {
62			TCanvas *can = new TCanvas("can","can");
63			can->SetLogy(1);
64
65			TH1F *FSR = PhotonSamples.Draw("FSR","genFSRNparticles" ,11,-0.5,10.5,"N(#gamma)", "events", BaseCut,mc,PlottingSetup::luminosity);
66			TH1F *all = PhotonSamples.Draw("all","genPhotonsNPhotons",11,-0.5,10.5,"N(#gamma)", "events", BaseCut,mc,PlottingSetup::luminosity);
67
68			FSR->SetLineColor(TColor::GetColor("#01DF01"));
69			all->SetLineColor(TColor::GetColor("#FF4000"));
70
71			TLegend *leg = new TLegend(0.7,0.7,0.89,0.89);
72			leg->SetFillColor(kWhite);
73			leg->AddEntry(FSR,"FSR photons","l");
74			leg->AddEntry(all,"all photons","l");
75
76			FSR->Draw("histo");
77			all->Draw("same,histo");
78			leg->SetHeader("DY MC:");
79			leg->Draw("same,histo");
80			DrawMCPrelim();
81
82			CompleteSave(can,"PhotonStudies/NPhotons");
83
84
85			TH1F *PhotonPt[10];
86			for(int i=0;i<10;i++) {
87			PhotonPt[i] = PhotonSamples.Draw("FSR","genPhotonsPt",25,0,100,"N(#gamma)", "events", FSRcut && TCut(("genFSRNparticles=="+any2string(i)).c_str()),mc,PlottingSetup::luminosity);
88			}
89
90			int DrawCounter=0;
91
92			TLegend *leg2 = make_legend();
93			for(int i=0;i<10;i++) {
94			cout << PhotonPt[i]->Integral() << endl;
95			if(PhotonPt[i]->Integral()>50) {
96			PhotonPt[i]->SetLineColor(DrawCounter+1);
97			leg2->AddEntry(PhotonPt[i],(any2string(i)+" FSR Photons").c_str(),"l");
98			if(DrawCounter) PhotonPt[i]->DrawNormalized("same,hist");
99			else PhotonPt[i]->DrawNormalized("hist");
100			DrawCounter++;
101			}
102			}
103
104			leg2->Draw();
105
106			leg2->SetHeader("DY MC (N_{jets} #geq 0):");
107			DrawMCPrelim();
108
109			CompleteSave(can,"PhotonStudies/PhotonPtvsNPhotons");
110
111			delete can;
112			}
113
114			void GenPhotons::MakeCharacterizationPlot(string variable, int nbins, float min, float max, string label, string saveas) {
115			TCanvas *can = new TCanvas("can","can");
116			can->SetLogy(1);
117
118			TH1F *dR_FSR = PhotonSamples.Draw("dR_FSR", variable, nbins, min, max, label, "events",BaseCut && FSRcut, mc,PlottingSetup::luminosity);
119			TH1F *dR_nFSR = PhotonSamples.Draw("dR_nFSR", variable, nbins, min, max, label, "events",BaseCut && NoFSRcut, mc,PlottingSetup::luminosity);
120
121			dR_FSR->SetLineColor(TColor::GetColor("#01DF01"));
122			dR_nFSR->SetLineColor(TColor::GetColor("#FF4000"));
123	buchmann	1.1
124	buchmann	1.2	TLegend *leg = new TLegend(0.7,0.7,0.89,0.89);
125			leg->SetFillColor(kWhite);
126			leg->AddEntry(dR_FSR,"FSR #gamma's","l");
127			leg->AddEntry(dR_nFSR,"non-FSR #gamma's","l");
128
129			if(dR_nFSR->GetMinimum()<dR_FSR->GetMinimum()) dR_FSR->SetMinimum(dR_nFSR->GetMinimum());
130			if(dR_nFSR->GetMaximum()>dR_FSR->GetMaximum()) dR_FSR->SetMaximum(dR_nFSR->GetMaximum());
131
132
133			dR_FSR->Draw("histo");
134			dR_nFSR->Draw("same,histo");
135			leg->SetHeader("DY MC:");
136			leg->Draw("same,histo");
137			DrawMCPrelim();
138
139			CompleteSave(can,"PhotonStudies/CharacterizingFSRPhotons/"+saveas);
140			delete dR_FSR;
141			delete dR_nFSR;
142			delete can;
143			}
144
145
146			void GenPhotons::CharacterizingFSRPhotons() {
147			MakeCharacterizationPlot("min(genPhotonsdR1,genPhotonsdR2)",100,0,5,"min(#Delta R (#gamma,l_{1}), #Delta R (#gamma,l_{2}))", "MinDR");
148			MakeCharacterizationPlot("genPhotonsPt",100,0,100,"p_{T}^{#gamma}", "PhotonPt");
149			MakeCharacterizationPlot("genPhotonsEta",100,-5,5,"{#eta}^{#gamma}", "PhotonEta");
150			MakeCharacterizationPlot("genPhotonsPhi",100,0,4,"{#phi}^{#gamma}", "PhotonPhi");
151			MakeCharacterizationPlot("genPhotonsM",100,0,2,"m^{#gamma}", "PhotonM");
152			}
153
154			void GenPhotons::CompareGenJZBDistributions() {
155			TCanvas *can = new TCanvas("can","can");
156			can->SetLogy(1);
157
158			TH1F *dR_FSR = PhotonSamples.Draw("dR_FSR", "genFSRjzbG", 30, -30, 30, "JZB^{gen}", "events",BaseCut && TCut("abs(genFSRjzbG)>0.1"), mc,PlottingSetup::luminosity);
159			TH1F *dR_nFSR = PhotonSamples.Draw("dR_nFSR", "genFSRjzb", 30, -30, 30, "JZB^{gen}", "events",BaseCut && TCut("abs(genFSRjzb)>0.1") , mc,PlottingSetup::luminosity);
160
161			dR_FSR->SetLineColor(TColor::GetColor("#01DF01"));
162			dR_nFSR->SetLineColor(TColor::GetColor("#FF4000"));
163
164			TLegend *leg = new TLegend(0.7,0.7,0.89,0.89);
165			leg->SetFillColor(kWhite);
166			leg->AddEntry(dR_FSR,"With FSR recovery","l");
167			leg->AddEntry(dR_nFSR,"Without FSR recovery","l");
168
169			if(dR_nFSR->GetMinimum()<dR_FSR->GetMinimum()) dR_FSR->SetMinimum(dR_nFSR->GetMinimum());
170			if(dR_nFSR->GetMaximum()>dR_FSR->GetMaximum()) dR_FSR->SetMaximum(dR_nFSR->GetMaximum());
171
172
173			dR_FSR->Draw("histo");
174			dR_nFSR->Draw("same,histo");
175			leg->SetHeader("DY MC:");
176			leg->Draw("same,histo");
177			DrawMCPrelim();
178
179			CompleteSave(can,"PhotonStudies/GenJZBDistributions");
180			delete dR_FSR;
181			delete dR_nFSR;
182			delete can;
183			}
184
185			void GenPhotons::dR2DPlots() {
186			TCanvas *can = new TCanvas("can","can",1000,500);
187			can->Divide(2,1);
188			can->cd(1);
189
190			vector<float> binX, binY;
191			for(int i=0;i<=31;i++) {
192			binX.push_back(i*0.1);
193			binY.push_back(i*0.1);
194			}
195
196			TH2F *dR_FSR = PhotonSamples.Draw("dR_FSR", "genPhotonsdR1:genPhotonsdR2", binX,binY, "#Delta R (#gamma^{FSR},l_{1})","#Delta R (#gamma^{FSR},l_{2})",BaseCut && FSRcut ,mc,PlottingSetup::luminosity, PhotonSamples.FindSample("o"),false);
197			TH2F *dR_nFSR = PhotonSamples.Draw("dR_nFSR", "genPhotonsdR1:genPhotonsdR2", binX,binY, "#Delta R (#gamma^{non-FSR},l_{1})","#Delta R (#gamma^{non-FSR},l_{2})",BaseCut && NoFSRcut ,mc,PlottingSetup::luminosity, PhotonSamples.FindSample("o"),false);
198
199			can->cd(1);
200			dR_FSR->Draw("COLZ");
201			DrawMCPrelim();
202
203			can->cd(2);
204			dR_nFSR->Draw("COLZ");
205			DrawMCPrelim();
206
207			CompleteSave(can,"PhotonStudies/dR_in_2D");
208			delete dR_FSR;
209			delete dR_nFSR;
210			delete can;
211			}
212
213
214			void GenPhotons::EfficiencyPurity() {
215			TCanvas *can = new TCanvas("can","can");
216			cout << "Computing efficiency and purity !" << endl;
217			vector<float> bindR, binPt;
218			for(int i=0;i<=100;i++) bindR.push_back(i*0.05);
219			binPt.push_back(0);binPt.push_back(10);binPt.push_back(50);binPt.push_back(100);binPt.push_back(1000);binPt.push_back(1010);
220			// binPt.push_back(0);binPt.push_back(1000);binPt.push_back(1001);
221
222
223	buchmann	1.3	TGraph *purity[binPt.size()];
224			TGraph *efficiency[binPt.size()];
225			TGraph *purity_n_efficiency[binPt.size()];
226			TH1F *FSR[binPt.size()];
227			TH1F *Any[binPt.size()];
228			for(int i=1;i<binPt.size();i++) {
229			//do this for each pt range
230			float ptlow=binPt[i-1];
231			float pthi =binPt[i];
232			stringstream ptcut;
233			ptcut << "genPhotonsPt>" << ptlow << " && genPhotonsPt<" << pthi << " && min(genPhotonsdR1,genPhotonsdR2)<5";
234
235			FSR[i-1] = PhotonSamples.Draw("FSR", "min(genPhotonsdR1,genPhotonsdR2)", bindR,"min DR","events",TCut(ptcut.str().c_str()) && BaseCut && FSRcut ,mc,PlottingSetup::luminosity);
236			Any[i-1] = PhotonSamples.Draw("Any", "min(genPhotonsdR1,genPhotonsdR2)", bindR,"min DR","events",TCut(ptcut.str().c_str()) && BaseCut ,mc,PlottingSetup::luminosity);
237
238			float NPhotons=0;
239			float NFSRPhotons=0;
240			float totalNFSRphotons=0;
241
242			for(int k=1;k<FSR[i-1]->GetNbinsX();k++) totalNFSRphotons+=FSR[i-1]->GetBinContent(k);
243			cout << i << " : " << totalNFSRphotons << " FSR photons in pt range " << ptlow << " , " << pthi << endl;
244			purity[i-1]=new TGraph(FSR[i-1]->GetNbinsX());
245			efficiency[i-1]=new TGraph(FSR[i-1]->GetNbinsX());
246			purity_n_efficiency[i-1]=new TGraph(FSR[i-1]->GetNbinsX());
247			for(int j=1;j<=FSR[i-1]->GetNbinsX();j++) {
248			NPhotons+=(int)Any[i-1]->GetBinContent(j);
249			NFSRPhotons+=FSR[i-1]->GetBinContent(j);
250			// cout << " ::::::::" << j << " : " << NPhotons << " : " << NFSRPhotons << " : " << totalNFSRphotons << endl;
251			if(NPhotons>0) purity[i-1]->SetPoint(j-1,FSR[i-1]->GetXaxis()->GetBinCenter(j)+FSR[i-1]->GetBinWidth(j),(float)NFSRPhotons/NPhotons);
252			else purity[i-1]->SetPoint(j-1,FSR[i-1]->GetXaxis()->GetBinCenter(j)+FSR[i-1]->GetBinWidth(j),-1);
253			if(totalNFSRphotons>0) efficiency[i-1]->SetPoint(j,FSR[i-1]->GetXaxis()->GetBinCenter(j)+FSR[i-1]->GetBinWidth(j),(float)NFSRPhotons/totalNFSRphotons);
254			else efficiency[i-1]->SetPoint(j-1,FSR[i-1]->GetXaxis()->GetBinCenter(j)+FSR[i-1]->GetBinWidth(j),-1);
255			if(totalNFSRphotons>0&&NPhotons>0) purity_n_efficiency[i-1]->SetPoint(j-1,(float)NFSRPhotons/totalNFSRphotons,(float)NFSRPhotons/NPhotons);
256	buchmann	1.2	}
257			cout << "Done for pt range" << endl;
258			}
259			cout << "Done with graphs" << endl;
260	buchmann	1.3	// cout << "Artificially aborting ... " << endl;return;
261			for(int i=0;i<binPt.size()-1;i++) {
262			cout << "About to draw purity for i=" << i << endl;
263			cout << "( this corresponds to the Pt range from " << binPt[i] << " to " << binPt[i+1] << ")" << endl;
264	buchmann	1.2	if(i==0) purity[i]->Draw("AC");
265			else purity[i]->Draw("C");
266			}
267			CompleteSave(can,"PhotonStudies/Purity");
268
269	buchmann	1.3	for(int i=0;i<binPt.size()-1;i++) {
270			cout << "About to draw efficiency for i=" << i << endl;
271	buchmann	1.2	if(i==0) efficiency[i]->Draw("AC");
272			else efficiency[i]->Draw("C");
273			}
274			CompleteSave(can,"PhotonStudies/Efficiency");
275
276	buchmann	1.3	for(int i=0;i<binPt.size()-1;i++) {
277	buchmann	1.2	if(i==0) purity_n_efficiency[i]->Draw("AC");
278			else purity_n_efficiency[i]->Draw("C");
279			}
280			CompleteSave(can,"PhotonStudies/Purity_n_efficiency");
281			}
282
283	buchmann	1.1
284			void GenPhotons::GenLevelStudies() {
285	buchmann	1.2	TCut essential_bkp = essentialcut;
286			essentialcut=TCut("");
287
288	buchmann	1.3	// CompareZmassToFinalStateLeptons();
289			// NumberOfPhotonsPerEvent();
290			// CharacterizingFSRPhotons();
291			// CompareGenJZBDistributions();
292			// dR2DPlots();
293	buchmann	1.2
294			EfficiencyPurity();
295	buchmann	1.1
296	buchmann	1.2	essentialcut=essential_bkp;
297	buchmann	1.1	}
298
299			void InitializePhotonSamples() {
300			float ZJetsCrossSection = 3532.8;
301			float LowZJetsCrossSection = 110500.0691.15; // LO xs * filter eff * k-factor (same as ZJets)
302
303	buchmann	1.2	PhotonSamples.AddSample("/shome/buchmann/ntuples/MC8tev/PhotonStudies/Skim3genJets_DYJetsToLL_M-50_TuneZ2Star_8TeV-madgraph-tarball-Summer12_DR53X-PU_S10_START53_V7A-v1.root","Z+Jets",-30346766,ZJetsCrossSection,false,false,1,kYellow);
304			PhotonSamples.AddSample("/shome/buchmann/ntuples/MC8tev/PhotonStudies/Skim3genJets_DYJetsToLL_M-10To50filter_8TeV-madgraph-Summer12_DR53X-PU_S10_START53_V7A-v1.root","Z+Jets",-6955271,LowZJetsCrossSection,false,false,1,kYellow);
305	buchmann	1.1
306			write_info(__FUNCTION__,"Photon samples have been initialized!");
307			PhotonSamples.ListSamples();
308
309	buchmann	1.2	GenPhotons::FSRcut = TCut("genPhotonsIsFSR");
310			GenPhotons::NoFSRcut = TCut("!genPhotonsIsFSR");
311			GenPhotons::BaseCut = TCut("genNjets>=3 && genPhotonsPt > 5 && NgenLeps>0");
312	buchmann	1.1	}
313
314			void FSRstudy() {
315			InitializePhotonSamples();
316
317			GenPhotons::GenLevelStudies();
318			RecoPhotons::RecoLevelStudies();
319			}