45 string cal_name[] = {
"iHCal",
"oHCal",
"EMCal" };
46 string cal_tag[] = {
"_in",
"_out",
"_emc" };
64 double dPhi = phi1 - phi2;
65 while (dPhi>=M_PI) { dPhi -= 2.*M_PI; }
66 while (dPhi<=-M_PI) { dPhi += 2.*M_PI; }
70 double dR(
double phi1,
double phi2,
double eta1,
double eta2){
72 double deta_sqrd = (eta1-eta2)*(eta1-eta2);
73 return sqrt(dphi_sqrd+deta_sqrd);
77 while (phi>=M_PI) { phi -= 2.*M_PI; }
78 while (phi<=-M_PI) { phi += 2.*M_PI; }
83 double min_eta = eta - 0.024 - 0.012;
84 double max_eta = eta + 0.024 + 0.012;
85 if ( calo_eta<min_eta || calo_eta>max_eta ) {
return false; }
86 double min_phi = phi - (M_PI/128.) - (M_PI/256.);
87 double max_phi = phi + (M_PI/128.) + (M_PI/256.);
88 if ( calo_phi<min_phi || calo_phi>max_phi ) {
return false; }
89 if ( calo_eta>=min_eta && calo_eta<=max_eta && calo_phi>=min_phi && calo_phi<=max_phi ) {
return true; }
93 TH1D *
hGeantPhi =
new TH1D(
"hGeantPhi",
";truth #phi",64,-M_PI,M_PI);
94 TH2D *
hEMCal3x3 =
new TH2D(
"hEMCal3x3",
"pion p_{T};E_{EMCal}^{3x3} about max-E tower [GeV]",30,0.,30.,120,0.,60.);
95 TH2D *
hEMCal3x3_FINE =
new TH2D(
"hEMCal3x3_FINE",
"pion p_{T};E_{EMCal}^{3x3} about max-E tower [GeV]",30,0.,30.,200.,0.,2.);
96 TH3D *
hE_inner_vs_outer =
new TH3D(
"hE_inner_vs_outer",
";pion p_{T};total E deposited in iHCal;total E deposited in oHCal",60,0.,30.,100, 0., 10.,120,0.,60.);
97 TH3D *
hDeltaPhi_fraction_HcalOverAll =
new TH3D(
"hDeltaPhi_fraction_HcalOverAll",
";pion p_{T};E_{HCals}/E_{all calos};oHCal #Delta#phi",60,0.,30.,100,0.,1.,160,-0.4,0.4);
98 TH3D *
hDeltaPhi_fraction_oHcalOverHcals =
new TH3D(
"hDeltaPhi_fraction_oHcalOverHcals",
";pion p_{T};E_{oHCal}/E_{HCals};oHCal #Delta#phi",60,0.,30.,100,0.,1.,160,-0.4,0.4);
108 hE_weighted_eta_phi[
i_cal] =
new TH3D(Form(
"hE_weighted_eta_phi_%s",
cal_name[
i_cal].c_str()),
";pion p_{T};calo #eta;calo #phi",60,0.,30.,
n_eta_bins[
i_cal],-
eta_max[i_cal],
eta_max[i_cal],
n_phi_bins[i_cal],-M_PI,M_PI);
109 h_eta_phi[
i_cal] =
new TH3D(Form(
"h_eta_phi_%s",
cal_name[i_cal].c_str()),
";pion p_{T};calo #eta;calo #phi",60,0.,30.,
n_eta_bins[i_cal],-
eta_max[i_cal],
eta_max[i_cal],
n_phi_bins[i_cal],-M_PI,M_PI);
111 hDeltaPhi_pt[
i_cal] =
new TH2D(Form(
"hDeltaPhi_pt_%s",
cal_name[i_cal].c_str()),
";pion p_{T};tower #phi - truth #phi",120,0.,60.,320,-4.,4.);
113 hDeltaPhi_fraction[
i_cal] =
new TH2D(Form(
"hDeltaPhi_fraction_%s",
cal_name[i_cal].c_str()),
";tower #phi - truth #phi;fraction of total calo E in layer",320,-4.,4.,100,0,1.);
114 hDeltaPhi_iPhi[
i_cal] =
new TH2D(Form(
"hDeltaPhi_iPhi_%s",
cal_name[i_cal].c_str()),
";tower #phi - truth #phi;tower iphi",320,-4.,4.,256,0.,256.);
118 hEnergy_fraction[
i_cal] =
new TH2D(Form(
"hEnergy_fraction_%s",
cal_name[i_cal].c_str()),
";pion E;fraction of pion E in calo",120,0.,60.,100,0.,10.);
119 hDeltaPhi[
i_cal] =
new TH1D(Form(
"hDeltaPhi_%s",
cal_name[i_cal].c_str()),
";tower #phi - truth #phi",128,-2.*M_PI,2.*M_PI);
125 TTree *
fChain = (TTree*)f->Get(
"T");
127 fChain->SetBranchAddress(
"event", &
event);
128 fChain->SetBranchAddress(
"eta", &
eta);
129 fChain->SetBranchAddress(
"phi", &
phi);
130 fChain->SetBranchAddress(
"e", &
e);
131 fChain->SetBranchAddress(
"pt", &
pt);
132 fChain->SetBranchAddress(
"vx", &
vx);
133 fChain->SetBranchAddress(
"vy", &
vy);
134 fChain->SetBranchAddress(
"vz", &
vz);
138 fChain->SetBranchAddress(Form(
"eta%s",
cal_tag[i_cal].c_str()), &
calo_eta[i_cal]);
139 fChain->SetBranchAddress(Form(
"phi%s",
cal_tag[i_cal].c_str()), &
calo_phi[i_cal]);
140 fChain->SetBranchAddress(Form(
"e%s",
cal_tag[i_cal].c_str()), &
calo_e[i_cal]);
141 fChain->SetBranchAddress(Form(
"ieta%s",
cal_tag[i_cal].c_str()), &
calo_ieta[i_cal]);
142 fChain->SetBranchAddress(Form(
"iphi%s",
cal_tag[i_cal].c_str()), &
calo_iphi[i_cal]);
145 Long64_t nentries = fChain->GetEntriesFast();
147 for (Long64_t jentry=0; jentry<nentries;jentry++) {
148 Long64_t ientry = fChain->GetTreeNumber();
149 if (ientry < 0)
break;
151 fChain->GetEntry(jentry);
159 float closestOuter_tow_dR = 9999.;
160 float closestOuter_tow_eta;
161 float closestOuter_tow_phi;
163 float energyIn3x3 = 0;
166 float max_tow_e = 0.;
171 for (
int itow=0; itow<
calo_e[
i_cal]->size(); ++itow) {
189 if (delta_R<closestOuter_tow_dR){
190 closestOuter_tow_dR = delta_R;
207 if (
i_cal==1) {max_out_phi=max_tow_phi;}
253 string remove_this =
".root";
254 size_t pos = outputroot.find(remove_this);
255 if (pos != string::npos) { outputroot.erase(pos, remove_this.length()); }
258 TFile *outFile =
new TFile(outFileName,
"RECREATE");
289 TCanvas * can =
new TCanvas(
"can" ,
"" ,700 ,500 );
291 const char us[] =
"_";
293 auto fa1 =
new TF1(
"fa1",
"0.",0,1);
295 for (
int i=0;
i<histos.size(); ++
i ) {
296 TH2D *hTemp = (TH2D*)histos[
i]->Project3D(
"ZY");
297 hTemp->SetLineColor(kRed);
298 hTemp->SetMarkerColor(kRed);
300 hTemp->ProfileX()->Draw(
"SAME");
301 can->SaveAs(Form(
"plots/ExploreTTrees/%s.pdf",histos[
i]->GetName()),
"PDF");
302 hTemp->ProfileX()->Draw(
"");
304 can->SaveAs(Form(
"plots/ExploreTTrees/%s_pfx.pdf",histos[
i]->GetName()),
"PDF");
305 for (
int j=0;
j<3; ++
j) {
307 double pthi = (
double)10.*(
j+1);
308 histos[
i]->GetXaxis()->SetRangeUser(ptlo,pthi);
309 hTemp = (TH2D*)histos[
i]->Project3D(
"ZY");
310 hTemp->SetLineColor(kRed);
311 hTemp->SetMarkerColor(kRed);
313 hTemp->ProfileX()->Draw(
"SAME");
314 can->SaveAs(Form(
"plots/ExploreTTrees/%s%s%i%s%i.pdf",histos[
i]->GetName(),us,(
int)ptlo,us,(
int)pthi),
"PDF");
315 hTemp->ProfileX()->Draw(
"");
317 can->SaveAs(Form(
"plots/ExploreTTrees/%s%s%i%s%i_pfx.pdf",histos[
i]->GetName(),us,(
int)ptlo,us,(
int)pthi),
"PDF");
319 for (
int j=1;
j<10; ++
j) {
321 double pthi = (
double)1.*(
j+1);
322 histos[
i]->GetXaxis()->SetRangeUser(ptlo,pthi);
323 hTemp = (TH2D*)histos[
i]->Project3D(
"ZY");
324 hTemp->SetLineColor(kRed);
325 hTemp->SetMarkerColor(kRed);
327 hTemp->ProfileX()->Draw(
"SAME");
328 can->SaveAs(Form(
"plots/ExploreTTrees/%s%s%i%s%i.pdf",histos[
i]->GetName(),us,(
int)ptlo,us,(
int)pthi),
"PDF");
329 hTemp->ProfileX()->Draw(
"");
331 can->SaveAs(Form(
"plots/ExploreTTrees/%s%s%i%s%i_pfx.pdf",histos[
i]->GetName(),us,(
int)ptlo,us,(
int)pthi),
"PDF");
333 histos[
i]->GetXaxis()->SetRangeUser(0.,30.);
339 can->SaveAs(
"plots/ExploreTTrees/hEMCal3x3.pdf",
"PDF");
342 can->SaveAs(
"plots/ExploreTTrees/hEMCal3x3_FINE.pdf",
"PDF");
345 TF1 *
f1 =
new TF1(
"f1",
"gaus",0.,0.6);
349 TH1D *hEMCal3x3_FINE_py[30];
351 for (
int i=1;
i<30; ++
i) {
355 hEMCal3x3_FINE_py[
i-1] = (TH1D*)
hEMCal3x3_FINE->ProjectionY(Form(
"hEMCal3x3_FINE_py__%i_%iGeV",ptlo,pthi),ptbin,ptbin);
357 hEMCal3x3_FINE_py[
i-1]->SetTitle(Form(
"%i-%i GeV/c pion",ptlo,pthi));
358 hEMCal3x3_FINE_py[
i-1]->SetStats(0);
359 hEMCal3x3_FINE_py[
i-1]->GetYaxis()->SetRangeUser(0.0,80.0);
360 hEMCal3x3_FINE_py[
i-1]->Draw(
"SAMEPLCPMC");
361 hEMCal3x3_FINE_py[
i-1]->Fit(f1,
"",
"",0.,0.6);
362 MIP_value[
i-1] = f1->GetMaximumX();
365 float mean_value = 0;
366 for (
int i=0;
i<29; ++
i) {
367 cout<<
i+1<<
"-"<<
i+2<<
"GeV/c \t"<<MIP_value[
i]<<endl;
368 mean_value += MIP_value[
i];
371 cout<<endl<<mean_value<<endl;
375 leg0 =
new TLegend(0.7, 0.3, 0.98, 0.98,NULL,
"brNDC");
376 leg0->SetBorderSize(0); leg0->SetLineColor(1); leg0->SetLineStyle(1);
377 leg0->SetLineWidth(1);
380 leg0->SetNColumns(2);
381 leg0->SetHeader(
"p_{T}^{pion}\t\t\tGaussian Mean",
"");
385 for (
int i=0;
i<29; ++
i) {
387 leg0->AddEntry(hEMCal3x3_FINE_py[
i],Form(
"%i-%i GeV/c",i+1,i+2),
"lpf");
388 leg0->AddEntry(
"",Form(
"%f",MIP_value[i]),
"");