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);
    
    
    TH2F *dR_FSR   = PhotonSamples.Draw("dR_FSR",  "genPhotonsPt:min(genPhotonsdR1,genPhotonsdR2)", bindR, binPt,"min DR","#gamma pT",BaseCut && FSRcut   ,mc,PlottingSetup::luminosity, PhotonSamples.FindSample("o"),false);
    TH2F *dR_Any  = PhotonSamples.Draw("dR_Any", "genPhotonsPt:min(genPhotonsdR1,genPhotonsdR2)",bindR,binPt, "min DR","#gamma pT",BaseCut  ,mc,PlottingSetup::luminosity, PhotonSamples.FindSample("o"),false);
    
    TGraph *purity[dR_FSR->GetNbinsY()+1];
    TGraph *efficiency[dR_FSR->GetNbinsY()+1];
    TGraph *purity_n_efficiency[dR_FSR->GetNbinsY()+1];
    for(int i=1;i<dR_FSR->GetNbinsY();i++) {
        //do this for each pt range
        int NPhotons=0;
        int NFSRPhotons=0;
        int totalNFSRphotons=0;
        for(int k=1;k<dR_FSR->GetNbinsX();k++) totalNFSRphotons+=dR_FSR->GetBinContent(k,i);
        cout << i << " : " << totalNFSRphotons << " FSR photons in pt range " << dR_FSR->GetYaxis()->GetBinLowEdge(i) << " , " << dR_FSR->GetYaxis()->GetBinLowEdge(i)+dR_FSR->GetYaxis()->GetBinWidth(i) << endl;
        purity[i-1]=new TGraph(dR_FSR->GetNbinsX());
        efficiency[i-1]=new TGraph(dR_FSR->GetNbinsX());
        purity_n_efficiency[i-1]=new TGraph(dR_FSR->GetNbinsX());
        for(int j=1;j<=dR_FSR->GetNbinsX();j++) {
            NPhotons+=dR_Any->GetBinContent(j,i);
            NFSRPhotons+=dR_FSR->GetBinContent(j,i);

//            cout << "For the pt range " << dR_FSR->GetYaxis()->GetBinLowEdge(i) << " , " << dR_FSR->GetYaxis()->GetBinLowEdge(i)+dR_FSR->GetYaxis()->GetBinWidth(i) << " we find a purity of " << (float)NFSRPhotons/NPhotons << " and an efficiency of " << (float)NFSRPhotons/totalNFSRphotons << endl;
            cout << "       ::::::::" << j << " : " << NPhotons << " : " << NFSRPhotons << " : " << totalNFSRphotons << endl;
/*            if(NPhotons>0) purity[i-1]->SetPoint(j-1,dR_FSR->GetXaxis()->GetBinCenter(j)+dR_FSR->GetBinWidth(j),(float)NFSRPhotons/NPhotons);
            else purity[i-1]->SetPoint(j-1,dR_FSR->GetXaxis()->GetBinCenter(j)+dR_FSR->GetBinWidth(j),-1);
            if(totalNFSRphotons>0) efficiency[i-1]->SetPoint(j,dR_FSR->GetXaxis()->GetBinCenter(j)+dR_FSR->GetBinWidth(j),(float)NFSRPhotons/totalNFSRphotons);
            else efficiency[i-1]->SetPoint(j-1,dR_FSR->GetXaxis()->GetBinCenter(j)+dR_FSR->GetBinWidth(j),-1);
         //   if(NPhotons>0 && totalNFSRphotons>0) purity_n_efficiency->SetPoint(j-1,(float)NFSRPhotons/totalNFSRphotons,(float)NFSRPhotons/NPhotons);
         //   else purity_n_efficiency->SetPoint(j-1,-1,-1);*/
        }
        cout << "Done for pt range" << endl;
    }
    cout << "Done with graphs" << endl;
    TFile *f = new TFile("out.root","RECREATE");
    dR_Any->Write();
    dR_FSR->Write();
    f->Close();
    cout << "Artificially aborting ... " << endl;return;
    for(int i=0;i<dR_FSR->GetNbinsY();i++) {
        if(i==0) purity[i]->Draw("AC");
        else purity[i]->Draw("C");
    }
    CompleteSave(can,"PhotonStudies/Purity");
    
    for(int i=0;i<dR_FSR->GetNbinsY();i++) {
        if(i==0) efficiency[i]->Draw("AC");
        else efficiency[i]->Draw("C");
    }
    CompleteSave(can,"PhotonStudies/Efficiency");
    
    for(int i=0;i<dR_FSR->GetNbinsY();i++) {
        if(i==0) purity_n_efficiency[i]->Draw("AC");
        else purity_n_efficiency[i]->Draw("C");
    }
    CompleteSave(can,"PhotonStudies/Purity_n_efficiency");
    
    dR_FSR->Draw("COLZ");
    CompleteSave(can,"PhotonStudies/dR_Pt_Diag_FSR");
    dR_Any->Draw("COLZ");
    CompleteSave(can,"PhotonStudies/dR_Pt_Diag_Any");
    
    
}


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.2
Committed:	Mon Oct 15 11:53:03 2012 UTC (12 years, 6 months ago) by buchmann
Content type:	text/plain
Branch:	MAIN
Changes since 1.1:	+262 -11 lines
Log Message:	Added a bunch of generator level FSR studies
#	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			TH2F *dR_FSR = PhotonSamples.Draw("dR_FSR", "genPhotonsPt:min(genPhotonsdR1,genPhotonsdR2)", bindR, binPt,"min DR","#gamma pT",BaseCut && FSRcut ,mc,PlottingSetup::luminosity, PhotonSamples.FindSample("o"),false);
224			TH2F *dR_Any = PhotonSamples.Draw("dR_Any", "genPhotonsPt:min(genPhotonsdR1,genPhotonsdR2)",bindR,binPt, "min DR","#gamma pT",BaseCut ,mc,PlottingSetup::luminosity, PhotonSamples.FindSample("o"),false);
225
226			TGraph *purity[dR_FSR->GetNbinsY()+1];
227			TGraph *efficiency[dR_FSR->GetNbinsY()+1];
228			TGraph *purity_n_efficiency[dR_FSR->GetNbinsY()+1];
229			for(int i=1;i<dR_FSR->GetNbinsY();i++) {
230			//do this for each pt range
231			int NPhotons=0;
232			int NFSRPhotons=0;
233			int totalNFSRphotons=0;
234			for(int k=1;k<dR_FSR->GetNbinsX();k++) totalNFSRphotons+=dR_FSR->GetBinContent(k,i);
235			cout << i << " : " << totalNFSRphotons << " FSR photons in pt range " << dR_FSR->GetYaxis()->GetBinLowEdge(i) << " , " << dR_FSR->GetYaxis()->GetBinLowEdge(i)+dR_FSR->GetYaxis()->GetBinWidth(i) << endl;
236			purity[i-1]=new TGraph(dR_FSR->GetNbinsX());
237			efficiency[i-1]=new TGraph(dR_FSR->GetNbinsX());
238			purity_n_efficiency[i-1]=new TGraph(dR_FSR->GetNbinsX());
239			for(int j=1;j<=dR_FSR->GetNbinsX();j++) {
240			NPhotons+=dR_Any->GetBinContent(j,i);
241			NFSRPhotons+=dR_FSR->GetBinContent(j,i);
242
243			// cout << "For the pt range " << dR_FSR->GetYaxis()->GetBinLowEdge(i) << " , " << dR_FSR->GetYaxis()->GetBinLowEdge(i)+dR_FSR->GetYaxis()->GetBinWidth(i) << " we find a purity of " << (float)NFSRPhotons/NPhotons << " and an efficiency of " << (float)NFSRPhotons/totalNFSRphotons << endl;
244			cout << " ::::::::" << j << " : " << NPhotons << " : " << NFSRPhotons << " : " << totalNFSRphotons << endl;
245			/* if(NPhotons>0) purity[i-1]->SetPoint(j-1,dR_FSR->GetXaxis()->GetBinCenter(j)+dR_FSR->GetBinWidth(j),(float)NFSRPhotons/NPhotons);
246			else purity[i-1]->SetPoint(j-1,dR_FSR->GetXaxis()->GetBinCenter(j)+dR_FSR->GetBinWidth(j),-1);
247			if(totalNFSRphotons>0) efficiency[i-1]->SetPoint(j,dR_FSR->GetXaxis()->GetBinCenter(j)+dR_FSR->GetBinWidth(j),(float)NFSRPhotons/totalNFSRphotons);
248			else efficiency[i-1]->SetPoint(j-1,dR_FSR->GetXaxis()->GetBinCenter(j)+dR_FSR->GetBinWidth(j),-1);
249			// if(NPhotons>0 && totalNFSRphotons>0) purity_n_efficiency->SetPoint(j-1,(float)NFSRPhotons/totalNFSRphotons,(float)NFSRPhotons/NPhotons);
250			// else purity_n_efficiency->SetPoint(j-1,-1,-1);*/
251			}
252			cout << "Done for pt range" << endl;
253			}
254			cout << "Done with graphs" << endl;
255			TFile *f = new TFile("out.root","RECREATE");
256			dR_Any->Write();
257			dR_FSR->Write();
258			f->Close();
259			cout << "Artificially aborting ... " << endl;return;
260			for(int i=0;i<dR_FSR->GetNbinsY();i++) {
261			if(i==0) purity[i]->Draw("AC");
262			else purity[i]->Draw("C");
263			}
264			CompleteSave(can,"PhotonStudies/Purity");
265
266			for(int i=0;i<dR_FSR->GetNbinsY();i++) {
267			if(i==0) efficiency[i]->Draw("AC");
268			else efficiency[i]->Draw("C");
269			}
270			CompleteSave(can,"PhotonStudies/Efficiency");
271
272			for(int i=0;i<dR_FSR->GetNbinsY();i++) {
273			if(i==0) purity_n_efficiency[i]->Draw("AC");
274			else purity_n_efficiency[i]->Draw("C");
275			}
276			CompleteSave(can,"PhotonStudies/Purity_n_efficiency");
277
278			dR_FSR->Draw("COLZ");
279			CompleteSave(can,"PhotonStudies/dR_Pt_Diag_FSR");
280			dR_Any->Draw("COLZ");
281			CompleteSave(can,"PhotonStudies/dR_Pt_Diag_Any");
282
283
284			}
285
286	buchmann	1.1
287			void GenPhotons::GenLevelStudies() {
288	buchmann	1.2	TCut essential_bkp = essentialcut;
289			essentialcut=TCut("");
290
291	buchmann	1.1	CompareZmassToFinalStateLeptons();
292	buchmann	1.2	NumberOfPhotonsPerEvent();
293			CharacterizingFSRPhotons();
294			CompareGenJZBDistributions();
295			dR2DPlots();
296
297			EfficiencyPurity();
298	buchmann	1.1
299	buchmann	1.2	essentialcut=essential_bkp;
300	buchmann	1.1	}
301
302			void InitializePhotonSamples() {
303			float ZJetsCrossSection = 3532.8;
304			float LowZJetsCrossSection = 110500.0691.15; // LO xs * filter eff * k-factor (same as ZJets)
305
306	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);
307			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);
308	buchmann	1.1
309			write_info(__FUNCTION__,"Photon samples have been initialized!");
310			PhotonSamples.ListSamples();
311
312	buchmann	1.2	GenPhotons::FSRcut = TCut("genPhotonsIsFSR");
313			GenPhotons::NoFSRcut = TCut("!genPhotonsIsFSR");
314			GenPhotons::BaseCut = TCut("genNjets>=3 && genPhotonsPt > 5 && NgenLeps>0");
315	buchmann	1.1	}
316
317			void FSRstudy() {
318			InitializePhotonSamples();
319
320			GenPhotons::GenLevelStudies();
321			RecoPhotons::RecoLevelStudies();
322			}