isoCluster/offline/analyzeClusterEtIso.C

0001 #include <sPhenixStyle.h>
0002 #include <sPhenixStyle.C>
0003
0004 /*
0005   This file analyzes the output of the primary
0006   isoCluster module, specifically the correlation
0007   between the cluster tranverse energy and the
0008   likelihood of clusters to be direct photons
0009   based on their Et values
0010
0011   Anthony Hodges, UIUC, February 15th, 2023
0012 */
0013
0014 void analyzeClusterEtIso(int pass)
0015 {
0016
0017   SetsPhenixStyle();
0018
0019   TFile *fin = new TFile(Form("../macro/condor/aggregation/gammaJet_pass%d.root",pass));
0020
0021   TTree *T = (TTree*)fin -> Get("T");
0022
0023   vector<float> *clusterPhi = 0;
0024   vector<float> *clusterEta = 0;
0025   vector<float> *clusterE = 0;
0026   vector<float> *clusterPt = 0;
0027   vector<float> *clusterEtIso = 0;
0028   vector<float> *clusterchisq = 0;
0029   vector<float> *photonPhi = 0;
0030   vector<float> *photonEta = 0;
0031   vector<float> *photonE = 0;
0032   vector<float> *photonPt = 0;
0033
0034   T -> SetBranchAddress("clusterPhi",&clusterPhi);
0035   T -> SetBranchAddress("clusterEta",&clusterEta);
0036   T -> SetBranchAddress("clusterE",&clusterE);
0037   T -> SetBranchAddress("clusterPt",&clusterPt);
0038   T -> SetBranchAddress("clusterEtIso",&clusterEtIso);
0039   T -> SetBranchAddress("photonPhi",&photonPhi);
0040   T -> SetBranchAddress("photonEta",&photonEta);
0041   T -> SetBranchAddress("photonE",&photonE);
0042   T -> SetBranchAddress("photonPt",&photonPt);
0043
0044   int nEvents = T -> GetEntries();
0045
0046   std::cout << "Analyzing " << nEvents << " entries" << std::endl;
0047
0048   TH2F *isoEtE = new TH2F("isoEtE","isoEtE",500,20,100,500,-10,10);
0049   isoEtE -> SetTitle(";E [GeV]; E_T [GeV]");
0050   TH2F *isoEtMatchFrac = new TH2F("isoEtMatchFrac","isoEtMatchFrac",500,-50,50,100,0,1);
0051   isoEtMatchFrac -> SetTitle("E_T [GeV]; Purity");
0052
0053   //eTIso calibration. Determine sigmalized values for Et based on clusterE
0054   float eBins[] = {25,30,35,40,45,50,55,60};
0055   int nBins = sizeof(eBins)/sizeof(eBins[0]);
0056   float fitStart = -60; float fitEnd = 60;
0057
0058   for(int i = 0; i < nEvents; i++)
0059     {
0060       T -> GetEntry(i);
0061
0062       for(int j = 0; j < clusterEtIso -> size(); j++)
0063     {
0064       isoEtE -> Fill(clusterE -> at(j), clusterEtIso -> at(j));
0065
0066     }
0067     }
0068
0069   TCanvas *c1 = new TCanvas();
0070   isoEtE -> Draw("colz");
0071
0072   TGraphErrors *sigma = new TGraphErrors();
0073   TGraphErrors *mu = new TGraphErrors();
0074
0075   for(int i = 0; i < nBins-1; i++)
0076     {
0077       TH1F *isoEtEProj = (TH1F*)isoEtE -> ProjectionY("projection",isoEtE -> GetXaxis() -> FindBin(eBins[i] + 0.001), isoEtE -> GetXaxis() -> FindBin(eBins[i+1] - 0.001));
0078
0079       TF1 *func = new TF1("func","gaus",fitStart, fitEnd);
0080       func -> SetLineColor(2);
0081       isoEtEProj -> Fit(func,"Q0","",fitStart,fitEnd);
0082       TCanvas *cTemp = new TCanvas();
0083
0084       isoEtEProj -> Fit(func,"Q","",func -> GetParameter(1) - func -> GetParameter(2), func -> GetParameter(1) + func -> GetParameter(2));
0085       sigma -> SetPoint(i, (eBins[i+1] + eBins[i])/2, func -> GetParameter(2));
0086       sigma -> SetPointError(i, 0, func -> GetParError(2));
0087
0088       mu ->  SetPoint(i, (eBins[i+1] + eBins[i])/2, func -> GetParameter(1));
0089       mu -> SetPointError(i, 0, func -> GetParError(1));
0090     }
0091
0092   TCanvas *c2 = new TCanvas();
0093   sigma -> SetTitle(";Cluster E [GeV]; Iso E_{T} #sigma [GeV]");
0094   sigma -> Draw("ap");
0095   TF1 *sigmaFit = new TF1("sigmaFit","pol1",eBins[0],eBins[nBins-1]);
0096   sigmaFit -> SetLineColor(2);
0097   sigma -> Fit(sigmaFit,"RQ","");
0098
0099   TCanvas *c3 = new TCanvas();
0100   mu -> SetTitle(";Cluster E [GeV];Iso E_{T} #mu [GeV]");
0101   mu -> Draw("ap");
0102   //TF1 *muFit = new TF1("muFit","pol8",eBins[0],eBins[nBins-1]);
0103   TF1 *muFit = new TF1("muFit","[0]*exp([1]*pow(x,[2]) + [3]) + [4]",(eBins[0] + eBins[1])/2,(eBins[nBins-2]+eBins[nBins-1])/2);
0104   muFit -> SetParameter(0, 0.00017);
0105
0106
0107   muFit -> SetLineColor(2);
0108   mu -> Fit(muFit,"RQ","");
0109   //Now we're going to go through the clusters again and assess
0110   //the correlation between eT and direct photon purity
0111   float sigmas[] = {0.5, 1 , 1.5, 2, 2.5, 3};
0112   const int nSigmas = sizeof(sigmas)/sizeof(sigmas[0]);
0113   int colors[] = {1,2,4,kGreen+2, kViolet,kCyan,kOrange+2,kMagenta+2,kAzure-2};
0114
0115   TGraphErrors *gPurity[nSigmas];
0116   TGraphErrors *gEfficiency[nSigmas];
0117
0118   TCanvas *c4 = new TCanvas();
0119   TCanvas *c5 = new TCanvas();
0120   TH1F *hdrMin = new TH1F("drMin","drMin",100,0,0.6);
0121   float dRMax = 0.12;
0122
0123   TH1F *hEResp = new TH1F("eResp","eResp",200,0,2);
0124   float respMin = 0.1;
0125   TLegend *leg = new TLegend(0.5,0.5,0.9,0.9);
0126   leg -> SetFillStyle(0);
0127   leg -> SetBorderSize(0);
0128
0129   for(int s = 0; s < nSigmas; s++)
0130     {
0131       gPurity[s] = new TGraphErrors();
0132       gPurity[s] -> SetMarkerColor(colors[s]);
0133
0134       gEfficiency[s] = new TGraphErrors();
0135       gEfficiency[s] -> SetMarkerColor(colors[s]);
0136       for(int e = 0; e < nBins; e++)
0137     {
0138
0139       int nMatch = 0;
0140       int nIsoCluster = 0;
0141       int nDirPhotons = 0;
0142       for(int i = 0; i < nEvents; i++)
0143         {
0144           T -> GetEntry(i);
0145
0146           if(photonE -> size() == 0) continue;
0147           nDirPhotons += photonE -> size();
0148           int matchIDCluster = -1;
0149           int matchIDPhoton = -1;
0150           int checkOnce = 0;
0151           for(int k = 0; k < clusterE -> size(); k++)
0152         {
0153           if(clusterE -> at(k) < eBins[e] || clusterE -> at(k) > eBins[e+1]) continue;
0154           if(abs(clusterEtIso -> at(k)) > muFit -> Eval(clusterE -> at(k)) + sigmas[s]*sigmaFit -> Eval(clusterE -> at(k))) continue;
0155           nIsoCluster++;
0156           float drMin = 1000;
0157
0158           for(int l = 0; l < photonE -> size(); l++)
0159             {
0160               if((photonE -> at(l) <  eBins[e] || photonE-> at(l) > eBins[e+1]) && !checkOnce)
0161             {
0162               checkOnce = 1;
0163               nDirPhotons++;
0164             }
0165               float dPhi = clusterPhi -> at(k) - photonPhi -> at(l);
0166               float dEta = clusterEta -> at(k) - photonEta -> at(l);
0167
0168               while(dPhi > M_PI) dPhi -= 2*M_PI;
0169               while(dPhi < -M_PI) dPhi += 2*M_PI;
0170
0171               float dr = sqrt(pow(dPhi,2) + pow(dEta,2));
0172
0173               if(dr < drMin)
0174             {
0175               drMin = dr;
0176               matchIDCluster = k;
0177               matchIDPhoton = l;
0178             }
0179             }//direct photon loop
0180
0181           hdrMin -> Fill(drMin);
0182
0183           if(!matchIDCluster || !matchIDPhoton) continue;
0184
0185           float response =  clusterE -> at(matchIDCluster)/photonE -> at(matchIDPhoton);
0186           hEResp -> Fill(response);
0187
0188           if(drMin < dRMax && abs(1 - response) < respMin) nMatch++;
0189
0190         }//cluster loop
0191
0192
0193         }//event loop
0194
0195       float purity = 0;
0196       if(nIsoCluster != 0) purity = (float)nMatch/(float)nIsoCluster;
0197       float efficiency = (float)nMatch/(float)nDirPhotons;
0198
0199
0200       // std::cout << Form("nMatch for sigma %g: %d with %d nIsoCluster",sigmas[s],nMatch,nIsoCluster) << std::endl;
0201       // std::cout << "filling with purity: " << purity << std::endl;
0202       gPurity[s] -> SetPoint(e, (eBins[e] + eBins[e+1])/2, purity);
0203       gEfficiency[s] -> SetPoint(e, (eBins[e] + eBins[e+1])/2, efficiency);
0204
0205
0206     }//energy loop
0207       leg -> AddEntry(gPurity[s], Form("E_{T} Iso #sigma < %g",sigmas[s]),"p");
0208       if(s == 0)
0209     {
0210       c4 -> cd();
0211       gPurity[s] -> Draw("ap");
0212       gPurity[s] -> SetTitle(";Cluster E [GeV]; Purity");
0213       gPurity[s] -> GetYaxis() -> SetRangeUser(0,0.02);
0214
0215       c5 -> cd();
0216       gEfficiency[s] -> Draw("ap");
0217       gEfficiency[s] -> SetTitle(";Cluster E [GeV]; Efficiency");
0218       gEfficiency[s] -> GetYaxis() -> SetRangeUser(0,0.02);
0219     }
0220       else
0221     {
0222       c4 -> cd();
0223       gPurity[s] -> Draw("samep");
0224       c5 -> cd();
0225       gEfficiency[s] -> Draw("samep");
0226     }
0227     }//sigma loop
0228   leg -> Draw("same");
0229
0230   TCanvas *cDrMin = new TCanvas();
0231   hdrMin -> Draw("epx0");
0232   hdrMin -> GetYaxis() -> SetRangeUser(1e-1, 10*hdrMin -> GetMaximum());
0233   TLine *ldrMax = new TLine(dRMax, 1e-1, dRMax, hdrMin -> GetMaximum());
0234   ldrMax -> SetLineStyle(7);
0235   ldrMax -> SetLineColor(2);
0236   ldrMax -> Draw("same");
0237   gPad -> SetLogy();
0238
0239   TCanvas *cResp = new TCanvas();
0240   hEResp -> Draw("epx0");
0241   TLine *lRespMin1 = new TLine(1-respMin, 0, 1-respMin, hEResp -> GetMaximum());
0242   lRespMin1 -> SetLineStyle(7);
0243   lRespMin1 -> SetLineColor(2);
0244   lRespMin1 -> Draw("same");
0245
0246   TLine *lRespMin2 = new TLine(1+respMin, 0, 1+respMin, hEResp -> GetMaximum());
0247   lRespMin2 -> SetLineStyle(7);
0248   lRespMin2 -> SetLineColor(2);
0249   lRespMin2 -> Draw("same");
0250
0251 }
0252
0253
0254
0255
0256
0257
0258
0259
0260
0261
0262