Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:12:14 

0001 // $Id: $
0002 
0003 /*!
0004  * \file Draw_PHG4DSTReader.C
0005  * \brief
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $Revision:   $
0008  * \date $Date: $
0009  */
0010 
0011 #include <cmath>
0012 #include <TFile.h>
0013 #include <TString.h>
0014 #include <TLine.h>
0015 #include <TTree.h>
0016 #include <cassert>
0017 #include "SaveCanvas.C"
0018 #include "SetOKStyle.C"
0019 using namespace std;
0020 
0021 TFile * _infile = NULL;
0022 TTree * T = NULL;
0023 TString cuts = "";
0024 
0025 void
0026 Draw_PHG4DSTReader( //
0027 //          const TString infile = "G4Hits_sPHENIX_e-_eta0_8GeV.root_DSTReader.root"
0028 //        const TString infile = "G4Hits_sPHENIX_e-_eta0_8GeV.root_Ana.root_DSTReader.root"//
0029 //    const TString infile =
0030 //        "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/single_particle/spacal2d/zerofield/G4Hits_sPHENIX_gamma_etaALL_50GeV.root_DSTReader.root" //
0031 //    const TString infile = "G4sPHENIXCells.root_DSTReader.root" //
0032 //    const TString infile = "G4sPHENIXCells_2DSpacal_100e10GeV.root_Ana.root_DSTReader.root"//
0033 //    const TString infile =
0034 //         "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_8GeV_CALICEBirk/Sim1096_eneg_DSTReader.root" //
0035 //    const TString infile =
0036 //        "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_8GeV/Sum_muonneg.lst_Process.root_DSTReader.root" //
0037 //    const TString infile = "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/Spacal2D_Leakage/Job_ALL_e-_DSTReader.root"//
0038 //        const TString infile = "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/Spacal2D_InsPHENIX/Sim1393ALL_eneg_DSTReader.root"//
0039     )
0040 {
0041   SetOKStyle();
0042   gStyle->SetOptStat(0);
0043   gStyle->SetOptFit(1111);
0044   TVirtualFitter::SetDefaultFitter("Minuit2");
0045   gSystem->Load("libg4eval.so");
0046 
0047   if (!_infile)
0048     {
0049       TString chian_str = infile;
0050       chian_str.ReplaceAll("ALL", "*");
0051 
0052       TChain * t = new TChain("T");
0053       const int n = t->Add(chian_str);
0054 
0055       cout << "Loaded " << n << " root files with " << chian_str << endl;
0056       assert(n>0);
0057 
0058       T = t;
0059 
0060       _infile = new TFile;
0061       _infile->SetName(infile);
0062     }
0063 
0064   assert(_infile);
0065 
0066   T->SetAlias("CEMC_Sample",
0067       "Sum$(G4HIT_CEMC.light_yield)/(Sum$(G4HIT_CEMC.edep) + Sum$(G4HIT_ABSORBER_CEMC.edep))");
0068   T->SetAlias("HCALIN_Sample",
0069       "Sum$(G4HIT_HCALIN.light_yield)/(Sum$(G4HIT_HCALIN.edep) + Sum$(G4HIT_ABSORBER_HCALIN.edep))");
0070   T->SetAlias("HCALOUT_Sample",
0071       "Sum$(G4HIT_HCALOUT.light_yield)/(Sum$(G4HIT_HCALOUT.edep) + Sum$(G4HIT_ABSORBER_HCALOUT.edep))");
0072   T->SetAlias("FHCAL_Sample",
0073       "Sum$(G4HIT_FHCAL.light_yield)/(Sum$(G4HIT_FHCAL.edep) + Sum$(G4HIT_ABSORBER_FHCAL.edep))");
0074 
0075   T->SetAlias("PHG4Particle0_e", "0+(PHG4Particle[0].fe)");
0076   T->SetAlias("PHG4Particle0_p",
0077       "1*sqrt(PHG4Particle[0].fpx**2+PHG4Particle[0].fpy**2+PHG4Particle[0].fpz**2)");
0078   T->SetAlias("PHG4Particle0_eta",
0079       "0.5*log((PHG4Particle0_p+PHG4Particle[0].fpz)/(PHG4Particle0_p-PHG4Particle[0].fpz))");
0080 
0081   T->SetAlias("PHG4Particle_p",
0082       "1*sqrt(PHG4Particle.fpx**2+PHG4Particle.fpy**2+PHG4Particle.fpz**2)");
0083   T->SetAlias("PHG4Particle_eta",
0084       "0.5*log((PHG4Particle_p+PHG4Particle.fpz)/(PHG4Particle_p-PHG4Particle.fpz))");
0085 
0086   T->SetAlias("Entrace_CEMC_x", "G4HIT_ABSORBER_CEMC[0].get_avg_x()+0");
0087   T->SetAlias("Entrace_CEMC_y", "G4HIT_ABSORBER_CEMC[0].get_avg_y()+0");
0088   T->SetAlias("Entrace_CEMC_z", "G4HIT_ABSORBER_CEMC[0].get_avg_z()+0");
0089   T->SetAlias("Entrace_CEMC_azimuthal",
0090       "atan2(Entrace_CEMC_y, Entrace_CEMC_x)*95.");
0091   T->SetAlias("Average_CEMC_z",
0092       "Sum$(G4HIT_ABSORBER_CEMC.get_avg_z()*G4HIT_ABSORBER_CEMC.edep)/Sum$(G4HIT_ABSORBER_CEMC.edep)");
0093   T->SetAlias("Average_BH_1_z",
0094       "Sum$(G4HIT_BH_1.get_avg_z())/Sum$(G4HIT_BH_1.hit_id>-1)");
0095 
0096   T->SetAlias("CEMC_E", "Sum$(G4HIT_CEMC.light_yield)+0");
0097   T->SetAlias("TOWER_CEMC_E",
0098       "Sum$(TOWER_CEMC.light_yield * (TOWER_CEMC.light_yield>32))");
0099   T->SetAlias("ABSORBER_CEMC_E", "Sum$(G4HIT_ABSORBER_CEMC.edep)+0");
0100   T->SetAlias("Total_CEMC_E",
0101       "Sum$(G4HIT_ABSORBER_CEMC.edep)+Sum$(G4HIT_CEMC.edep)");
0102   T->SetAlias("HCALIN_E", "Sum$(G4HIT_HCALIN.light_yield)+0");
0103   T->SetAlias("HCALOUT_E", "Sum$(G4HIT_HCALOUT.light_yield)+0");
0104   T->SetAlias("Total_HCALIN_E",
0105       "Sum$(G4HIT_ABSORBER_HCALIN.edep)+Sum$(G4HIT_HCALIN.edep)");
0106   T->SetAlias("Total_HCALOUT_E",
0107       "Sum$(G4HIT_ABSORBER_HCALOUT.edep)+Sum$(G4HIT_HCALOUT.edep)");
0108   T->SetAlias("BH_1_E", "Sum$(G4HIT_BH_1.edep)+0");
0109   T->SetAlias("EMCScale", "1*1");
0110 
0111   T->SetAlias("PHG4Particle0_pT",
0112       "1*sqrt(PHG4Particle[0].fpx**2+PHG4Particle[0].fpy**2)");
0113   T->SetAlias("PHG4Particle0_theta",
0114       "1*atan2(PHG4Particle0_pT,PHG4Particle[0].fpz)");
0115   T->SetAlias("PHG4Particle0_phi",
0116       "1*atan2(PHG4Particle[0].fpy, PHG4Particle[0].fpx)");
0117 
0118   T->SetAlias("PHG4Particle0_xhx",
0119       "cos(PHG4Particle0_theta)*cos(PHG4Particle0_phi)");
0120   T->SetAlias("PHG4Particle0_xhy",
0121       "cos(PHG4Particle0_theta)*sin(PHG4Particle0_phi)");
0122   T->SetAlias("PHG4Particle0_xhz", "-sin(PHG4Particle0_theta)*1");
0123   T->SetAlias("PHG4Particle0_yhx", "1*cos(PHG4Particle0_phi + pi/2)");
0124   T->SetAlias("PHG4Particle0_yhy", "1*sin(PHG4Particle0_phi + pi/2)");
0125   T->SetAlias("PHG4Particle0_yhz", "0*1");
0126   T->SetAlias("PHG4Particle0_zhz", "cos(PHG4Particle0_theta)*1");
0127   T->SetAlias("PHG4Particle0_zhy",
0128       "sin(PHG4Particle0_theta)*cos(PHG4Particle0_phi)");
0129   T->SetAlias("PHG4Particle0_zhz",
0130       "sin(PHG4Particle0_theta)*sin(PHG4Particle0_phi)");
0131 
0132   T->SetAlias("Field", "0*0");
0133   TString field_cor = "";
0134 //  T->SetAlias("Field", "1*1");
0135 //  TString field_cor = "_field_cor";
0136   T->SetAlias("CEMC_yp_cor",
0137       "(-19.984 + 5.64862 * PHG4Particle0_pT -0.531359* PHG4Particle0_pT * PHG4Particle0_pT) * (-Field) ");
0138   T->SetAlias("ABSORBER_CEMC_yp_cor", "CEMC_yp_cor + 0");
0139   T->SetAlias("HCALIN_yp_cor",
0140       "(-25.6969 + 6.24025 * PHG4Particle0_pT -0.472037 * PHG4Particle0_pT * PHG4Particle0_pT) * (-Field) ");
0141   T->SetAlias("HCALOUT_yp_cor",
0142       "(-40.6969 + 6.24025 * PHG4Particle0_pT -0.472037 * PHG4Particle0_pT * PHG4Particle0_pT) * (-Field) ");
0143 
0144   T->SetAlias("CEMC_xp",
0145       "G4HIT_CEMC.get_avg_x() * PHG4Particle0_xhx + G4HIT_CEMC.get_avg_y() * PHG4Particle0_xhy + G4HIT_CEMC.get_avg_z() * PHG4Particle0_xhz");
0146   T->SetAlias("CEMC_yp",
0147       "CEMC_yp_cor + G4HIT_CEMC.get_avg_x() * PHG4Particle0_yhx + G4HIT_CEMC.get_avg_y() * PHG4Particle0_yhy + G4HIT_CEMC.get_avg_z() * PHG4Particle0_yhz");
0148   T->SetAlias("CEMC_zp",
0149       "G4HIT_CEMC.get_avg_x() * PHG4Particle0_zhx + G4HIT_CEMC.get_avg_y() * PHG4Particle0_zhy + G4HIT_CEMC.get_avg_z() * PHG4Particle0_zhz");
0150   T->SetAlias("CEMC_R", "sqrt(CEMC_xp**2 + CEMC_yp**2)");
0151 
0152   T->SetAlias("ABSORBER_CEMC_xp",
0153       "G4HIT_ABSORBER_CEMC.get_avg_x() * PHG4Particle0_xhx + G4HIT_ABSORBER_CEMC.get_avg_y() * PHG4Particle0_xhy + G4HIT_ABSORBER_CEMC.get_avg_z() * PHG4Particle0_xhz");
0154   T->SetAlias("ABSORBER_CEMC_yp",
0155       "ABSORBER_CEMC_yp_cor + G4HIT_ABSORBER_CEMC.get_avg_x() * PHG4Particle0_yhx + G4HIT_ABSORBER_CEMC.get_avg_y() * PHG4Particle0_yhy + G4HIT_ABSORBER_CEMC.get_avg_z() * PHG4Particle0_yhz");
0156   T->SetAlias("ABSORBER_CEMC_zp",
0157       "G4HIT_ABSORBER_CEMC.get_avg_x() * PHG4Particle0_zhx + G4HIT_ABSORBER_CEMC.get_avg_y() * PHG4Particle0_zhy + G4HIT_ABSORBER_CEMC.get_avg_z() * PHG4Particle0_zhz");
0158   T->SetAlias("ABSORBER_CEMC_R",
0159       "sqrt(ABSORBER_CEMC_xp**2 + ABSORBER_CEMC_yp**2)");
0160 
0161   T->SetAlias("HCALIN_xp",
0162       "G4HIT_HCALIN.get_avg_x() * PHG4Particle0_xhx + G4HIT_HCALIN.get_avg_y() * PHG4Particle0_xhy + G4HIT_HCALIN.get_avg_z() * PHG4Particle0_xhz");
0163   T->SetAlias("HCALIN_yp",
0164       "HCALIN_yp_cor + G4HIT_HCALIN.get_avg_x() * PHG4Particle0_yhx + G4HIT_HCALIN.get_avg_y() * PHG4Particle0_yhy + G4HIT_HCALIN.get_avg_z() * PHG4Particle0_yhz");
0165   T->SetAlias("HCALIN_zp",
0166       "G4HIT_HCALIN.get_avg_x() * PHG4Particle0_zhx + G4HIT_HCALIN.get_avg_y() * PHG4Particle0_zhy + G4HIT_HCALIN.get_avg_z() * PHG4Particle0_zhz");
0167   T->SetAlias("HCALIN_R", "sqrt(HCALIN_xp**2 + HCALIN_yp**2)");
0168 
0169   T->SetAlias("HCALOUT_xp",
0170       "G4HIT_HCALOUT.get_avg_x() * PHG4Particle0_xhx + G4HIT_HCALOUT.get_avg_y() * PHG4Particle0_xhy + G4HIT_HCALOUT.get_avg_z() * PHG4Particle0_xhz");
0171   T->SetAlias("HCALOUT_yp",
0172       "HCALOUT_yp_cor + G4HIT_HCALOUT.get_avg_x() * PHG4Particle0_yhx + G4HIT_HCALOUT.get_avg_y() * PHG4Particle0_yhy + G4HIT_HCALOUT.get_avg_z() * PHG4Particle0_yhz");
0173   T->SetAlias("HCALOUT_zp",
0174       "G4HIT_HCALOUT.get_avg_x() * PHG4Particle0_zhx + G4HIT_HCALOUT.get_avg_y() * PHG4Particle0_zhy + G4HIT_HCALOUT.get_avg_z() * PHG4Particle0_zhz");
0175   T->SetAlias("HCALOUT_R", "sqrt(HCALOUT_xp**2 + HCALOUT_yp**2)");
0176 
0177   T->SetAlias("TOWER_CEMC_eta_mean",
0178       "Sum$(TOWER_CEMC.bineta*TOWER_CEMC.get_energy())/Sum$(TOWER_CEMC.get_energy())");
0179   T->SetAlias("G4HIT_CEMC_eta_mean",
0180       "Sum$(-log(tan(0.5*atan2(sqrt(G4HIT_CEMC.get_avg_x()*G4HIT_CEMC.get_avg_x() + G4HIT_CEMC.get_avg_y()*G4HIT_CEMC.get_avg_y()) , G4HIT_CEMC.get_avg_z())))  * G4HIT_CEMC.edep)/Sum$(G4HIT_CEMC.edep)");
0181 
0182   T->SetAlias("TOWER_CEMC_phi_mean",
0183       "Sum$(TOWER_CEMC.binphi*TOWER_CEMC.get_energy())/Sum$(TOWER_CEMC.get_energy())");
0184   T->SetAlias("G4HIT_CEMC_phi_mean",
0185       "Sum$(atan2(G4HIT_CEMC.get_avg_y(),G4HIT_CEMC.get_avg_x()) * G4HIT_CEMC.edep)/Sum$(G4HIT_CEMC.edep)");
0186 
0187   const TCut event_sel = "1*1";
0188   cuts = "_all_event";
0189 
0190   cout << "Build event selection of " << (const char *) event_sel << endl;
0191 
0192   T->Draw(">>EventList", event_sel);
0193   TEventList * elist = gDirectory->GetObjectChecked("EventList", "TEventList");
0194   cout << elist->GetN() << " / " << T->GetEntriesFast() << " events selected"
0195       << endl;
0196 
0197   T->SetEventList(elist);
0198 ////
0199 //  DrawDist(infile);
0200 //  DrawLeakage(infile);
0201 //  DrawLeakage_Wide(infile);
0202 //  DrawLeakage_Phi(infile);
0203   Sampling(infile);
0204 //  DrawTower_Raw_E(infile);
0205 
0206 //  DrawCalibratedE(infile, 0.02206);
0207 //  DrawCalibratedE(infile, 0.02206*9.71762904052642762e-01);
0208 //
0209 //    DrawCalibratedE(infile, 0.02206*7.71244e+00/7.86120e+00);
0210 
0211 //  DrawCalibratedE_Tower(infile, 500.);
0212 //  DrawCalibratedE_Tower(infile, 500.*9.71762904052642762e-01);
0213 //  DrawCalibratedE_Tower(infile, 500.*7.71244e+00/7.86120e+00);
0214 
0215 //  DrawLeakage_LY(infile);
0216 //
0217 //  DrawTowerIDCheck(infile);
0218 }
0219 
0220 void
0221 DrawCalibratedE(TString infile, const double SF = 0.021)
0222 {
0223 
0224   TCanvas *c1 = new TCanvas("DrawCalibratedE" + cuts, "DrawCalibratedE" + cuts,
0225       1800, 900);
0226   c1->Divide(2, 1);
0227   int idx = 1;
0228   TPad * p;
0229 
0230   p = (TPad *) c1->cd(idx++);
0231   c1->Update();
0232   T->Draw("CEMC_E>>hCEMC_E(1000,0,.5)");
0233 
0234   p = (TPad *) c1->cd(idx++);
0235   c1->Update();
0236   T->Draw(Form("CEMC_E/%f>>hCEMC_E_SF(1000,0,15)", SF));
0237 
0238   SaveCanvas(c1,
0239       TString(_infile->GetName()) + TString("_DrawJet_")
0240           + TString(c1->GetName()), kFALSE);
0241 }
0242 
0243 void
0244 DrawCalibratedE_Tower(TString infile, const double SF = 0.021)
0245 {
0246 
0247   TCanvas *c1 = new TCanvas("DrawCalibratedE" + cuts, "DrawCalibratedE" + cuts,
0248       1800, 900);
0249   c1->Divide(2, 1);
0250   int idx = 1;
0251   TPad * p;
0252 
0253   p = (TPad *) c1->cd(idx++);
0254   c1->Update();
0255   T->Draw("TOWER_CEMC_E>>hCEMC_E(1000,0,8000)");
0256 
0257   p = (TPad *) c1->cd(idx++);
0258   c1->Update();
0259   T->Draw(Form("TOWER_CEMC_E/%f>>hCEMC_E_SF(1000,0,15)", SF));
0260 
0261   SaveCanvas(c1,
0262       TString(_infile->GetName()) + TString("_DrawJet_")
0263           + TString(c1->GetName()), kFALSE);
0264 }
0265 
0266 void
0267 DrawTower_Raw_E(TString infile)
0268 {
0269   SetOKStyle();
0270   gStyle->SetOptStat(0);
0271   gStyle->SetOptFit(0);
0272   TVirtualFitter::SetDefaultFitter("Minuit2");
0273   gSystem->Load("libg4eval.so");
0274 
0275   TCanvas *c1 = new TCanvas("DrawTower_Raw_E" + cuts, "DrawTower_Raw_E" + cuts,
0276       1100, 900);
0277   c1->Divide(1, 1);
0278   int idx = 1;
0279   TPad * p;
0280 
0281   p = (TPad *) c1->cd(idx++);
0282   c1->Update();
0283   p->SetLogy();
0284   p->SetGridx(0);
0285   p->SetGridy(0);
0286 
0287   T->Draw("TOWER_RAW_CEMC.energy>>he(100,0,25000)");
0288 
0289   he->SetTitle(";Photo-electron count per tower; A. U.");
0290   he->SetLineWidth(4);
0291   he->SetLineColor(kBlue + 1);
0292 
0293   SaveCanvas(c1,
0294       TString(_infile->GetName()) + "_Draw_PHG4DSTReader_"
0295           + TString(c1->GetName()), true);
0296 }
0297 
0298 void
0299 DrawCalibratedE_PlotTestBeam(
0300     const TString base_name =
0301         "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_Data")
0302 {
0303   SetOKStyle();
0304   gStyle->SetOptStat(0);
0305   gStyle->SetOptFit(0);
0306   TVirtualFitter::SetDefaultFitter("Minuit2");
0307   gSystem->Load("libg4eval.so");
0308 
0309   const TString config = "nocut";
0310 
0311   const int N = 8;
0312   const double beam_res = 2.7e-2;
0313 
0314   const double beam_E[8] =
0315     { 0.995, 2.026, 3.0, 4.12, 6.0, 8.0, 10.0, 12.0 };
0316   double res[100] =
0317     { 0 };
0318   double res_sub[100] =
0319     { 0 };
0320   double res_err[100] =
0321     { 0 };
0322 
0323   TCanvas *c1 = new TCanvas(config, config);
0324   c1->Print(base_name + "/input_" + config + ".pdf[");
0325 
0326   for (int i = 0; i < N; i++)
0327     {
0328       TFile *f = new TFile(
0329           base_name + Form("/%.0fGeV_", beam_E[i]) + config + ".root");
0330       assert(f);
0331 
0332       TH1D * ECalSumElectron = (TH1D *) f->GetObjectChecked("ECalSumElectron",
0333           "TH1D");
0334       TH1D * ECalSumNonElectron = (TH1D *) f->GetObjectChecked(
0335           "ECalSumNonElectron", "TH1D");
0336       assert(ECalSumElectron);
0337       assert(ECalSumNonElectron);
0338 //      ECalSumElectron->SetStats(false);
0339 //      ECalSumNonElectron->SetStats(false);
0340       TF1 * fgaus =
0341       ECalSumElectron->GetListOfFunctions()->At(0);
0342       assert(fgaus);
0343 
0344       ECalSumElectron->Draw();
0345       c1->Print(base_name + "/input_" + config + ".pdf");
0346 
0347       res[i] = fgaus->GetParameter(2)/fgaus->GetParameter(1);
0348       res_sub[i] = sqrt(res[i]*res[i] - beam_res*beam_res);
0349       res_err[i] = fgaus->GetParError(1)/fgaus->GetParameter(1);
0350     }
0351   c1->Print(base_name + "/input_" + config + ".pdf]");
0352 
0353   TGraphErrors * gamma_eta0 = new TGraphErrors(N, beam_E, res_sub, 0, res_err);
0354   gamma_eta0->Print();
0355   TGraphErrors * gamma_eta9 = new TGraphErrors(N, beam_E, res, 0, res_err);
0356   gamma_eta0->Print();
0357 
0358   TCanvas *c1 = new TCanvas("DrawCalibratedE_PlotTestBeam", "DrawCalibratedE_PlotTestBeam",
0359       1100, 900);
0360   c1->Divide(1, 1);
0361   int idx = 1;
0362   TPad * p;
0363 
0364   p = (TPad *) c1->cd(idx++);
0365   c1->Update();
0366 
0367   p->SetGridx(0);
0368   p->SetGridy(0);
0369 
0370 
0371   p->DrawFrame(0, 0, 14, 20e-2,
0372       ";Beam Energy (GeV);Relative energy resolution, #DeltaE/E");
0373 
0374   TLegend * lg = new TLegend(2.*14./35., 9.6e-2, 17.*14./35., 19.6e-2, NULL, "br");
0375   TLegend * lg2 = new TLegend(4, 9e-2, 33.*14./35., 19e-2, NULL, "br");
0376 
0377   TF1 * f_calo = new TF1("f_calo_gamma_eta0", "sqrt([0]*[0]+[1]*[1]/x)/100",
0378       0.5, 60);
0379   TF1 * f_calo_l = new TF1("f_calo_l_gamma_eta0", "([0]+[1]/sqrt(x))/100", 0.5,
0380       60);
0381   gamma_eta0->Fit(f_calo, "RM0");
0382   f_calo->Print();
0383   gamma_eta0->Fit(f_calo_l, "RM0");
0384   f_calo_l->Print();
0385 
0386   TF1 * f_caloR = new TF1("f_calo_gamma_eta0R", "sqrt([0]*[0]+[1]*[1]/x)/100",
0387       1.5, 60);
0388   TF1 * f_calo_lR = new TF1("f_calo_l_gamma_eta0R", "([0]+[1]/sqrt(x))/100", 1.5,
0389       60);
0390   gamma_eta0->Fit(f_caloR, "RM0");
0391   f_caloR->Print();
0392   gamma_eta0->Fit(f_calo_lR, "RM0");
0393   f_calo_lR->Print();
0394 
0395 
0396   gamma_eta0->SetLineColor(kBlack);
0397   gamma_eta0->SetMarkerColor(kBlack);
0398   gamma_eta0->SetLineWidth(2);
0399   gamma_eta0->SetMarkerStyle(kFullSquare);
0400   gamma_eta0->SetMarkerSize(2);
0401 
0402   gamma_eta9->SetLineColor(kBlack);
0403   gamma_eta9->SetMarkerColor(kBlack);
0404   gamma_eta9->SetLineWidth(2);
0405   gamma_eta9->SetMarkerStyle(kOpenSquare);
0406   gamma_eta9->SetMarkerSize(2);
0407 
0408   f_calo->SetLineColor(kRed + 1);
0409   f_calo->SetLineWidth(2);
0410   f_calo_l->SetLineColor(kRed + 1);
0411   f_calo_l->SetLineWidth(2);
0412   f_calo_l->SetLineStyle(kDashed);
0413 
0414   f_caloR->SetLineColor(kBlue + 3);
0415   f_caloR->SetLineWidth(2);
0416   f_calo_lR->SetLineColor(kBlue + 3);
0417   f_calo_lR->SetLineWidth(2);
0418   f_calo_lR->SetLineStyle(kDashed);
0419 
0420   f_calo->Draw("same");
0421   f_calo_l->Draw("same");
0422   f_caloR->Draw("same");
0423   f_calo_lR->Draw("same");
0424   gamma_eta0->Draw("p");
0425   gamma_eta9->Draw("p");
0426 
0427   lg2->AddEntry(gamma_eta9,
0428       Form("Electrons Data, #eta = 0.3-0.4", abs(f_calo->GetParameter(0)),
0429           f_calo->GetParameter(1)), "p");
0430   lg2->AddEntry(gamma_eta0,
0431       Form("Electrons Data - 2.7%% Beam #DeltaE", f_calo->GetParameter(0),
0432           f_calo->GetParameter(1)), "p");
0433   lg2->AddEntry(f_calo,
0434       Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",abs( f_calo->GetParameter(0)),
0435           f_calo->GetParameter(1)), "l");
0436   TLegendEntry * entry = lg2->AddEntry(f_calo_l,
0437       Form("#DeltaE/E = %.1f%%  + %.1f%%/#sqrt{E}", abs(f_calo_l->GetParameter(0)),
0438           f_calo_l->GetParameter(1)), "l");
0439   entry->SetTextColor(kGray + 1);
0440   lg2->AddEntry(f_caloR,
0441       Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}, E #geq 2 GeV", abs(f_caloR->GetParameter(0)),
0442           f_caloR->GetParameter(1)), "l");
0443   TLegendEntry * entry = lg2->AddEntry(f_calo_lR,
0444       Form("#DeltaE/E = %.1f%%  + %.1f%%/#sqrt{E}, E #geq 2 GeV",abs( f_calo_lR->GetParameter(0)),
0445           f_calo_lR->GetParameter(1)), "l");
0446   entry->SetTextColor(kGray + 1);
0447 
0448 //  lg->Draw();
0449   lg2->Draw();
0450 
0451   SaveCanvas(c1, base_name + "/" + TString(c1->GetName()), true);
0452 }
0453 
0454 void
0455 //DrawCalibratedE_Plot(TString ID = "4.12GeV", TString config = "")
0456 //DrawCalibratedE_Plot(TString ID = "8GeV", TString config = "")
0457 DrawCalibratedE_Plot(TString ID = "12GeV",TString config = "")
0458 //DrawCalibratedE_Plot(TString ID = "8GeV",TString config = "_Range1um")
0459 //DrawCalibratedE_Plot(TString ID = "8GeV",TString config = "_CALICEBirk")
0460 //DrawCalibratedE_Plot(TString ID = "4GeV",TString config = "")
0461 //DrawCalibratedE_Plot(TString ID = "4GeV",TString config = "_Range1um")
0462 //DrawCalibratedE_Plot(TString ID = "4GeV",TString config = "_CALICEBirk")
0463 {
0464   SetOKStyle();
0465   gStyle->SetOptStat(0);
0466   gStyle->SetOptFit(0);
0467   TVirtualFitter::SetDefaultFitter("Minuit2");
0468   gSystem->Load("libg4eval.so");
0469 
0470 //  /direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_8GeV//SimALL_gamma_DSTReader.root_DrawJet_DrawCalibratedE_all_event.root
0471 
0472 //  TFile *f = new TFile("/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_"+ID+config+"//SimALL_eneg_DSTReader.root_DrawJet_DrawCalibratedE_all_event.root");
0473 //  assert(f);
0474 //  TH1F * h_eneg = (TH1F *)f->GetObjectChecked("hCEMC_E_SF","TH1F");
0475 //  assert(h_eneg);
0476 //  TFile *f = new TFile("/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_"+ID+config+"//SimALL_pineg_DSTReader.root_DrawJet_DrawCalibratedE_all_event.root");
0477 //  assert(f);
0478 //  TH1F * h_pineg = (TH1F *)f->GetObjectChecked("hCEMC_E_SF","TH1F");
0479 //  assert(h_pineg);
0480 //  TFile *f = new TFile("/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_"+ID+config+"//SimALL_kaonneg_DSTReader.root_DrawJet_DrawCalibratedE_all_event.root");
0481 //  assert(f);
0482 //  TH1F * h_kaonneg = (TH1F *)f->GetObjectChecked("hCEMC_E_SF","TH1F");
0483 //  assert(h_kaonneg);
0484 
0485   TFile *f =
0486       new TFile(
0487           "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_"
0488               + ID + config
0489               + "//Sum_eneg.lst_Process.root_DSTReader.root_DrawJet_DrawCalibratedE_all_event.root");
0490   assert(f);
0491   TH1F * h_eneg = (TH1F *) f->GetObjectChecked("hCEMC_E_SF", "TH1F");
0492   assert(h_eneg);
0493   TFile *f =
0494       new TFile(
0495           "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_"
0496               + ID + config
0497               + "//Sum_pineg.lst_Process.root_DSTReader.root_DrawJet_DrawCalibratedE_all_event.root");
0498   assert(f);
0499   TH1F * h_pineg = (TH1F *) f->GetObjectChecked("hCEMC_E_SF", "TH1F");
0500   assert(h_pineg);
0501   TFile *f =
0502       new TFile(
0503           "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_"
0504               + ID + config
0505               + "//Sum_kaonneg.lst_Process.root_DSTReader.root_DrawJet_DrawCalibratedE_all_event.root");
0506   assert(f);
0507   TH1F * h_kaonneg = (TH1F *) f->GetObjectChecked("hCEMC_E_SF", "TH1F");
0508   assert(h_kaonneg);
0509   TFile *f =
0510       new TFile(
0511           "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_"
0512               + ID + config
0513               + "//Sum_muonneg.lst_Process.root_DSTReader.root_DrawJet_DrawCalibratedE_all_event.root");
0514   assert(f);
0515   TH1F * h_muonneg = (TH1F *) f->GetObjectChecked("hCEMC_E_SF", "TH1F");
0516   assert(h_muonneg);
0517 
0518   TH1F * h_eneg_scale_hadron_data_tail = (TH1F *) h_eneg->Clone(
0519       "h_eneg_scale_hadron_data_tail");
0520 
0521   TFile *f =
0522       new TFile(
0523           "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_8GeV/beam_data/"
0524               + ID + ".root");
0525   assert(f);
0526 
0527   TH1D * ECalSumElectron = (TH1D *) f->GetObjectChecked("ECalSumElectron",
0528       "TH1D");
0529   TH1D * ECalSumNonElectron = (TH1D *) f->GetObjectChecked("ECalSumNonElectron",
0530       "TH1D");
0531   assert(ECalSumElectron);
0532   assert(ECalSumNonElectron);
0533   ECalSumElectron->SetStats(false);
0534   ECalSumNonElectron->SetStats(false);
0535   ECalSumElectron->GetListOfFunctions()->RemoveAt(0);
0536 
0537   const double sim_E = 7.86131e+00;
0538   const double data_E = 1.09505e+03;
0539 
0540   ECalSumElectron->GetXaxis()->Set(ECalSumElectron->GetXaxis()->GetNbins(),
0541       ECalSumElectron->GetXaxis()->GetXmin() * sim_E / data_E,
0542       ECalSumElectron->GetXaxis()->GetXmax() * sim_E / data_E //
0543           );
0544   ECalSumNonElectron->GetXaxis()->Set(
0545       ECalSumNonElectron->GetXaxis()->GetNbins(),
0546       ECalSumNonElectron->GetXaxis()->GetXmin() * sim_E / data_E,
0547       ECalSumNonElectron->GetXaxis()->GetXmax() * sim_E / data_E //
0548           );
0549 
0550   const double energy = h_eneg->GetMean();
0551 
0552   double eneg_scale = ECalSumElectron->Integral(
0553       ECalSumElectron->GetXaxis()->FindBin(energy * .7),
0554       ECalSumElectron->GetXaxis()->FindBin(energy * 1.3))
0555       / h_eneg->Integral(h_eneg->GetXaxis()->FindBin(energy * .7),
0556           h_eneg->GetXaxis()->FindBin(energy * 1.3));
0557 
0558   const double hadron_norm_max_range = 0.7;
0559 
0560   double muonneg_scale = 0.1 * //10% muon in hadron sample
0561       ECalSumNonElectron->Integral(ECalSumNonElectron->GetXaxis()->FindBin(.01),
0562           ECalSumNonElectron->GetXaxis()->FindBin(energy * 1.3))
0563       / h_muonneg->Integral(h_muonneg->GetXaxis()->FindBin(.01),
0564           h_muonneg->GetXaxis()->FindBin(energy * 1.3));
0565 
0566   const double muon_ratio_in_MIP = h_muonneg->Integral(
0567       h_muonneg->GetXaxis()->FindBin(.01),
0568       h_muonneg->GetXaxis()->FindBin(hadron_norm_max_range))
0569       / h_muonneg->Integral(h_muonneg->GetXaxis()->FindBin(.01),
0570           h_muonneg->GetXaxis()->FindBin(energy * 1.3));
0571   const double muon_count_in_MIP = muon_ratio_in_MIP //
0572   * 0.1 * //10% muon in hadron sample
0573       ECalSumNonElectron->Integral(ECalSumNonElectron->GetXaxis()->FindBin(.01),
0574           ECalSumNonElectron->GetXaxis()->FindBin(energy * 1.3));
0575 
0576   cout << "DrawCalibratedE_Plot: muon_ratio_in_MIP = " << muon_ratio_in_MIP
0577       << ", muon_count_in_MIP = " << muon_count_in_MIP << endl;
0578 
0579   double pineg_scale = (ECalSumNonElectron->Integral(
0580       ECalSumNonElectron->GetXaxis()->FindBin(.01),
0581       ECalSumNonElectron->GetXaxis()->FindBin(hadron_norm_max_range))
0582       - muon_count_in_MIP)
0583       / h_pineg->Integral(h_pineg->GetXaxis()->FindBin(.01),
0584           h_pineg->GetXaxis()->FindBin(hadron_norm_max_range));
0585 
0586   double kaonneg_scale = (ECalSumNonElectron->Integral(
0587       ECalSumNonElectron->GetXaxis()->FindBin(.01),
0588       ECalSumNonElectron->GetXaxis()->FindBin(hadron_norm_max_range))
0589       - muon_count_in_MIP)
0590       / h_kaonneg->Integral(h_kaonneg->GetXaxis()->FindBin(.01),
0591           h_kaonneg->GetXaxis()->FindBin(hadron_norm_max_range));
0592 
0593   double eneg_scale_hadron_data_tail = ECalSumNonElectron->Integral(
0594       ECalSumNonElectron->GetXaxis()->FindBin(energy * 0.95),
0595       ECalSumNonElectron->GetXaxis()->FindBin(energy * 1.15))
0596       / h_eneg_scale_hadron_data_tail->Integral(
0597           h_eneg_scale_hadron_data_tail->GetXaxis()->FindBin(energy * 0.95),
0598           h_eneg_scale_hadron_data_tail->GetXaxis()->FindBin(energy * 1.15));
0599 
0600   h_eneg->Rebin(5);
0601   h_pineg->Rebin(5);
0602   h_kaonneg->Rebin(5);
0603   h_muonneg->Rebin(5);
0604   h_eneg_scale_hadron_data_tail->Rebin(5);
0605 
0606 //  ECalSumElectron->Rebin();
0607   ECalSumNonElectron->Rebin(40);
0608 
0609   ECalSumElectron->Sumw2();
0610   ECalSumNonElectron->Sumw2();
0611 
0612   ECalSumElectron->SetMarkerStyle(kFullCircle);
0613   ECalSumElectron->SetMarkerSize(1);
0614   ECalSumNonElectron->SetMarkerStyle(kFullCircle);
0615   ECalSumNonElectron->SetMarkerSize(1);
0616 
0617   eneg_scale *= ECalSumElectron->GetBinWidth(1) / h_eneg->GetBinWidth(1);
0618   pineg_scale *= ECalSumNonElectron->GetBinWidth(1) / h_pineg->GetBinWidth(1);
0619   kaonneg_scale *= ECalSumNonElectron->GetBinWidth(1)
0620       / h_kaonneg->GetBinWidth(1);
0621   muonneg_scale *= ECalSumNonElectron->GetBinWidth(1)
0622       / h_muonneg->GetBinWidth(1);
0623   eneg_scale_hadron_data_tail *= ECalSumNonElectron->GetBinWidth(1)
0624       / h_eneg_scale_hadron_data_tail->GetBinWidth(1);
0625 
0626   //Resolution for electron sample
0627 
0628   TF1 * f_gaus_sim = new TF1("f_gaus_sim", "gaus",
0629       h_eneg->GetMean() - 2 * h_eneg->GetRMS(),
0630       h_eneg->GetMean() + 2 * h_eneg->GetRMS());
0631   f_gaus_sim->SetParameters(1, h_eneg->GetMean(), h_eneg->GetRMS());
0632   new TCanvas;
0633   h_eneg->Fit(f_gaus_sim, "RM0");
0634   h_eneg->DrawClone();
0635   f_gaus_sim->Draw("same");
0636 
0637   TF1 * f_gaus_data = new TF1("f_gaus_data", "gaus",
0638       h_eneg->GetMean() - 2 * h_eneg->GetRMS(),
0639       h_eneg->GetMean() + 2 * h_eneg->GetRMS());
0640   f_gaus_data->SetParameters(1, h_eneg->GetMean(), h_eneg->GetRMS());
0641   new TCanvas;
0642   ECalSumElectron->Fit(f_gaus_data, "RM0");
0643   ECalSumElectron->DrawClone();
0644   f_gaus_data->Draw("same");
0645 
0646   h_eneg->Scale(eneg_scale);
0647   h_pineg->Scale(pineg_scale);
0648   h_kaonneg->Scale(kaonneg_scale);
0649   h_muonneg->Scale(muonneg_scale);
0650   h_eneg_scale_hadron_data_tail->Scale(eneg_scale_hadron_data_tail);
0651 
0652 //  h_pineg -> Add(h_eneg_scale_hadron_data_tail);
0653 //  h_kaonneg -> Add(h_eneg_scale_hadron_data_tail);
0654 
0655   TCanvas *c1 = new TCanvas("DrawCalibratedE_Plot" + cuts,
0656       "DrawCalibratedE_Plot" + cuts, 1800, 700);
0657   c1->Divide(2, 1);
0658   int idx = 1;
0659   TPad * p;
0660 
0661   p = (TPad *) c1->cd(idx++);
0662   c1->Update();
0663   p->SetLogy();
0664   p->SetGridx(0);
0665   p->SetGridy(0);
0666 
0667   h_eneg->SetLineWidth(2);
0668 
0669   ECalSumElectron->GetXaxis()->SetRangeUser(0, energy * 1.5);
0670   ECalSumElectron->Draw();
0671   h_eneg->SetLineWidth(3);
0672   h_eneg->SetLineColor(kGreen + 3);
0673   h_eneg->SetFillColor(kGreen - 7);
0674   h_eneg->SetFillStyle(1001);
0675   h_eneg->Draw("same");
0676   ECalSumElectron->Draw("same");
0677 
0678   ECalSumElectron->SetTitle(";Calibrated Energy (GeV);Count/bin");
0679 
0680   TLegend * lg = new TLegend(0.13, 0.65, 0.55, 0.85);
0681   lg->AddEntry(ECalSumElectron,
0682       Form("#splitline{2014 eRD1 beam test (e-cut)}{#DeltaE/E = %.1f%%}",
0683           f_gaus_data->GetParameter(2) / f_gaus_data->GetParameter(1) * 100),
0684       "pe");
0685   lg->AddEntry(h_eneg,
0686       Form("#splitline{e^{-} sim using sPHENIX Geant4,}{#DeltaE/E = %.1f%%}",
0687           f_gaus_sim->GetParameter(2) / f_gaus_sim->GetParameter(1) * 100),
0688       "lf");
0689   lg->Draw();
0690 
0691   p = (TPad *) c1->cd(idx++);
0692   c1->Update();
0693   p->SetLogy();
0694   p->SetGridx(0);
0695   p->SetGridy(0);
0696 
0697   h_pineg->SetLineColor(kRed);
0698   h_pineg->SetLineWidth(3);
0699   h_kaonneg->SetLineColor(kBlue);
0700   h_kaonneg->SetLineWidth(3);
0701   h_muonneg->SetLineColor(kBlack);
0702   h_muonneg->SetFillColor(kGray);
0703   h_muonneg->SetLineWidth(3);
0704   h_muonneg->SetFillStyle(1001);
0705   h_eneg_scale_hadron_data_tail->SetLineColor(kGreen + 3);
0706   h_eneg_scale_hadron_data_tail->SetLineWidth(2);
0707   h_eneg_scale_hadron_data_tail->SetFillColor(kGreen - 7);
0708   h_eneg_scale_hadron_data_tail->SetFillStyle(1001);
0709 
0710   ECalSumNonElectron->GetXaxis()->SetRangeUser(0, energy * 1.5);
0711 
0712   ECalSumNonElectron->Draw();
0713   h_muonneg->Draw("same");
0714 //  h_eneg_scale_hadron_data_tail->Draw("same");  ////////////// Guess of electron contamination
0715   h_pineg->Draw("same");
0716   h_kaonneg->Draw("same");
0717   ECalSumNonElectron->Draw("same");
0718 
0719   ECalSumNonElectron->SetTitle(";Calibrated Energy (GeV);Count/bin");
0720 
0721   TLegend * lg = new TLegend(0.35, 0.65, 0.85, 0.9);
0722   lg->AddEntry(ECalSumNonElectron, Form("2014 eRD1 beam test (#bar{e-cut})"),
0723       "pe");
0724   lg->AddEntry(h_pineg, Form("#pi^{-} sim using sPHENIX Geant4"), "l");
0725   lg->AddEntry(h_kaonneg, Form("K^{-} sim using sPHENIX Geant4"), "l");
0726 //  lg->AddEntry(h_eneg_scale_hadron_data_tail, Form("e^{-} sim using sPHENIX Geant4"),////////////// Guess of electron contamination
0727 //      "lf");////////////// Guess of electron contamination
0728   lg->AddEntry(h_muonneg, Form("#mu^{-} (~10%% data) sim using sPHENIX Geant4"), "lf");
0729   lg->Draw();
0730 
0731   SaveCanvas(c1,
0732       TString(
0733           "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_"
0734               + ID + config + "/DrawCalibratedE_Plot_DrawJet_")
0735           + TString(c1->GetName()), true);
0736 }
0737 
0738 void
0739 DrawTowerIDCheck(TString infile)
0740 {
0741 
0742   TCanvas *c1 = new TCanvas("DrawTowerIDCheck" + cuts,
0743       "DrawTowerIDCheck" + cuts, 1800, 900);
0744   c1->Divide(2, 1);
0745   int idx = 1;
0746   TPad * p;
0747 
0748   p = (TPad *) c1->cd(idx++);
0749   c1->Update();
0750 
0751   T->Draw("TOWER_CEMC_eta_mean:G4HIT_CEMC_eta_mean", "", "*");
0752 
0753   p = (TPad *) c1->cd(idx++);
0754   c1->Update();
0755 
0756   T->Draw("TOWER_CEMC_phi_mean:G4HIT_CEMC_phi_mean", "", "*");
0757 
0758   SaveCanvas(c1,
0759       TString(_infile->GetName()) + TString("_DrawJet_")
0760           + TString(c1->GetName()), kFALSE);
0761 }
0762 
0763 void
0764 DrawDist(TString infile)
0765 {
0766 
0767   TCanvas *c1 = new TCanvas("DrawDist" + cuts, "DrawDist" + cuts, 1800, 900);
0768   c1->Divide(2, 1);
0769   int idx = 1;
0770   TPad * p;
0771 
0772   p = (TPad *) c1->cd(idx++);
0773   c1->Update();
0774 
0775   T->Draw(
0776       "G4HIT_CEMC.get_avg_y():G4HIT_CEMC.get_avg_x()>>h_EMC_Azim(600,-300,300,600,-300,300)",
0777       "G4HIT_CEMC.edep", "goff");
0778 //  T->Draw(
0779 //      "G4HIT_HCALIN.get_avg_y():G4HIT_HCALIN.get_avg_x()>>h_HCALIN_Azim(600,-300,300,600,-300,300)",
0780 //      "G4HIT_HCALIN.edep", "goff");
0781 //  T->Draw(
0782 //      "G4HIT_HCALOUT.get_avg_y():G4HIT_HCALOUT.get_avg_x()>>h_HCALOUT_Azim(600,-300,300,600,-300,300)",
0783 //      "G4HIT_HCALOUT.edep", "goff");
0784 //  TGraph * g_EMC =
0785 //      (TGraph *) gPad->GetListOfPrimitives()->FindObject("Graph")->Clone(
0786 //          "g_EMC");
0787 //  assert(g_EMC);
0788 //  g_EMC->SetName("g_EMC");
0789 //  g_EMC->SetMarkerColor(kRed);
0790 //  g_EMC->SetMarkerStyle(kDot);
0791 
0792   p->DrawFrame(-300, -300, 300, 300,
0793       "Azimuthal distribution of hits;X (cm);Y (cm)");
0794   h_EMC_Azim->Draw("colz same");
0795 //  h_HCALIN_Azim->Draw("colz same");
0796 //  h_HCALOUT_Azim->Draw("colz same");
0797 
0798   TEllipse * cBaBar = new TEllipse(0, 0, 140);
0799 
0800   cBaBar->SetFillStyle(0);
0801   cBaBar->SetLineWidth(3);
0802   cBaBar->SetLineColor(kMagenta);
0803   cBaBar->Draw();
0804   TEllipse * cBaBar = new TEllipse(0, 0, 170);
0805 
0806   cBaBar->SetFillStyle(0);
0807   cBaBar->SetLineWidth(3);
0808   cBaBar->SetLineColor(kMagenta);
0809   cBaBar->Draw();
0810 
0811   p = (TPad *) c1->cd(idx++);
0812   c1->Update();
0813 
0814   T->Draw(
0815       "sqrt(G4HIT_CEMC.get_avg_x()**2 + G4HIT_CEMC.get_avg_y()**2):G4HIT_CEMC.get_avg_z()>>h_CEMC_Rz(600,-300,300,300,0,300)",
0816       "G4HIT_CEMC.edep", "goff");
0817 //  T->Draw(
0818 //      "sqrt(G4HIT_HCALIN.get_avg_x()**2 + G4HIT_HCALIN.get_avg_y()**2):G4HIT_HCALIN.get_avg_z()>>h_HCALIN_Rz(600,-300,300,300,0,300)",
0819 //      "G4HIT_HCALIN.edep", "goff");
0820 //  T->Draw(
0821 //      "sqrt(G4HIT_HCALOUT.get_avg_x()**2 + G4HIT_HCALOUT.get_avg_y()**2):G4HIT_HCALOUT.get_avg_z()>>h_HCALOUT_Rz(600,-300,300,300,0,300)",
0822 //      "G4HIT_HCALOUT.edep", "goff");
0823 
0824 //  TGraph * g_EMC =
0825 //      (TGraph *) gPad->GetListOfPrimitives()->FindObject("Graph")->Clone(
0826 //          "g_EMC");
0827 //  assert(g_EMC);
0828 //  g_EMC->SetName("g_EMC");
0829 //  g_EMC->SetMarkerColor(kRed);
0830 //  g_EMC->SetMarkerStyle(kDot);
0831 
0832   p->DrawFrame(-300, 00, 300, 300,
0833       "Azimuthal distribution of hits;Z (cm);R (cm)");
0834   h_CEMC_Rz->Draw("colz same");
0835 //  h_HCALIN_Rz->Draw("colz same");
0836 //  h_HCALOUT_Rz->Draw("colz same");
0837 
0838   TLine * lBaBar = new TLine(-174.5, 140, 174.5, 140);
0839 
0840   lBaBar->SetLineWidth(3);
0841   lBaBar->SetLineColor(kMagenta);
0842   lBaBar->Draw();
0843   TLine * lBaBar = new TLine(-174.5, 170, 174.5, 170);
0844 
0845   lBaBar->SetLineWidth(3);
0846   lBaBar->SetLineColor(kMagenta);
0847   lBaBar->Draw();
0848 
0849   SaveCanvas(c1,
0850       TString(_infile->GetName()) + TString("_DrawJet_")
0851           + TString(c1->GetName()), kFALSE);
0852 }
0853 
0854 void
0855 DrawLeakage(TString infile)
0856 {
0857 
0858   TCanvas *c1 = new TCanvas("DrawLeakage" + cuts, "DrawLeakage" + cuts, 1800,
0859       900);
0860   c1->Divide(2, 3);
0861   int idx = 1;
0862   TPad * p;
0863 
0864   p = (TPad *) c1->cd(idx++);
0865   c1->Update();
0866   T->Draw(
0867       "Average_CEMC_z:PHG4Particle0_eta>>hAverage_CEMC_z_PHG4Particle0_eta(120,0.8,0.9,120,75,125)",
0868       "", "colz");
0869 
0870   p = (TPad *) c1->cd(idx++);
0871   c1->Update();
0872   T->Draw(
0873       "Entrace_CEMC_z:PHG4Particle0_eta>>hEntrace_CEMC_z_PHG4Particle0_eta(120,0.8,0.9,120,75,125)",
0874       "", "colz");
0875 
0876   p = (TPad *) c1->cd(idx++);
0877   c1->Update();
0878 
0879   T->Draw(
0880       "Total_HCALIN_E/PHG4Particle0_e:Average_CEMC_z>>hLeakage_Average_CEMC_z(120,75,125,100,0,0.25)",
0881       "", "colz");
0882 
0883   p = (TPad *) c1->cd(idx++);
0884   c1->Update();
0885 
0886   T->Draw(
0887       "Total_HCALIN_E/PHG4Particle0_e:PHG4Particle0_eta>>hLeakage_PHG4Particle0_eta(120,0.8,0.9,100,0,0.25)",
0888       "", "colz");
0889 
0890   p = (TPad *) c1->cd(idx++);
0891   c1->Update();
0892   T->Draw(
0893       "Total_CEMC_E/PHG4Particle0_e:Average_CEMC_z>>hTotal_CEMC_E_Entrance_Z_Entrace_CEMC_z(120,75,125,100,0,1.2)",
0894       "", "colz");
0895 
0896   p = (TPad *) c1->cd(idx++);
0897   c1->Update();
0898   T->Draw(
0899       "Total_CEMC_E/PHG4Particle0_e:PHG4Particle0_eta>>hTotal_CEMC_E_PHG4Particle0_eta(120,0.8,0.9,100,0,1.2)",
0900       "", "colz");
0901 
0902   SaveCanvas(c1,
0903       TString(_infile->GetName()) + TString("_DrawJet_")
0904           + TString(c1->GetName()), kFALSE);
0905 }
0906 void
0907 DrawLeakage_LY(TString infile)
0908 {
0909 
0910   TCanvas *c1 = new TCanvas("DrawLeakage_LY" + cuts, "DrawLeakage_LY" + cuts,
0911       1800, 900);
0912   c1->Divide(2, 3);
0913   int idx = 1;
0914   TPad * p;
0915 
0916   p = (TPad *) c1->cd(idx++);
0917   c1->Update();
0918   T->Draw(
0919       "Sum$(G4HIT_CEMC.light_yield)/PHG4Particle0_e:Average_CEMC_z>>hTotal_CEMC_E_Entrance_Z_Entrace_CEMC_z(120,75,125,100,0,.05)",
0920       "", "colz");
0921 
0922   p = (TPad *) c1->cd(idx++);
0923   c1->Update();
0924   T->Draw(
0925       "Sum$(G4HIT_CEMC.light_yield)/PHG4Particle0_e:PHG4Particle0_eta>>hTotal_CEMC_E_PHG4Particle0_eta(120,0.8,0.9,100,0,.05)",
0926       "", "colz");
0927 
0928   SaveCanvas(c1,
0929       TString(_infile->GetName()) + TString("_DrawJet_")
0930           + TString(c1->GetName()), kFALSE);
0931 }
0932 
0933 void
0934 DrawLeakage_Wide(TString infile)
0935 {
0936 
0937   TCanvas *c1 = new TCanvas("DrawLeakage_Wide" + cuts,
0938       "DrawLeakage_Wide" + cuts, 1800, 900);
0939   c1->Divide(2, 3);
0940   int idx = 1;
0941   TPad * p;
0942 
0943   p = (TPad *) c1->cd(idx++);
0944   c1->Update();
0945   T->Draw(
0946       "Average_CEMC_z:PHG4Particle0_eta>>hAverage_CEMC_z_PHG4Particle0_eta(300,0,1.1,300,0,150)",
0947       "", "colz");
0948 
0949   p = (TPad *) c1->cd(idx++);
0950   c1->Update();
0951   T->Draw(
0952       "Entrace_CEMC_z:PHG4Particle0_eta>>hEntrace_CEMC_z_PHG4Particle0_eta(300,0,1.1,300,0,150)",
0953       "", "colz");
0954 
0955   p = (TPad *) c1->cd(idx++);
0956   c1->Update();
0957 
0958   T->Draw(
0959       "Total_HCALIN_E/PHG4Particle0_e:Average_CEMC_z>>hLeakage_Average_CEMC_z(300,0,150,100,0,0.25)",
0960       "", "colz");
0961 
0962   p = (TPad *) c1->cd(idx++);
0963   c1->Update();
0964 
0965   T->Draw(
0966       "Total_HCALIN_E/PHG4Particle0_e:PHG4Particle0_eta>>hLeakage_PHG4Particle0_eta(300,0,1.1,100,0,0.25)",
0967       "", "colz");
0968 
0969   p = (TPad *) c1->cd(idx++);
0970   c1->Update();
0971   T->Draw(
0972       "Total_CEMC_E/PHG4Particle0_e:Average_CEMC_z>>hTotal_CEMC_E_Entrance_Z_Entrace_CEMC_z(300,0,150,100,0,1.2)",
0973       "", "colz");
0974 
0975   p = (TPad *) c1->cd(idx++);
0976   c1->Update();
0977   T->Draw(
0978       "Total_CEMC_E/PHG4Particle0_e:PHG4Particle0_eta>>hTotal_CEMC_E_PHG4Particle0_eta(300,0,1.1,100,0,1.2)",
0979       "", "colz");
0980 
0981   SaveCanvas(c1,
0982       TString(_infile->GetName()) + TString("_DrawJet_")
0983           + TString(c1->GetName()), kFALSE);
0984 }
0985 
0986 void
0987 DrawLeakage_Phi(TString infile)
0988 {
0989 
0990   TCanvas *c1 = new TCanvas("DrawLeakage_Phi" + cuts, "DrawLeakage_Phi" + cuts,
0991       1800, 900);
0992   c1->Divide(1, 2);
0993   int idx = 1;
0994   TPad * p;
0995 
0996   p = (TPad *) c1->cd(idx++);
0997   c1->Update();
0998 
0999   T->Draw(
1000       "Total_HCALIN_E/PHG4Particle0_e:PHG4Particle0_phi>>hLeakage_PHG4Particle0_e(320,-3.14159,3.14159,100,0,0.25)",
1001       "", "colz");
1002 
1003   p = (TPad *) c1->cd(idx++);
1004   c1->Update();
1005   T->Draw(
1006       "Total_CEMC_E/PHG4Particle0_e:PHG4Particle0_phi>>hTotal_CEMC_E_PHG4Particle0_e(320,-3.14159,3.14159,100,0,1.2)",
1007       "", "colz");
1008 
1009   SaveCanvas(c1,
1010       TString(_infile->GetName()) + TString("_DrawJet_")
1011           + TString(c1->GetName()), kFALSE);
1012 }
1013 
1014 void
1015 DrawLeakage_Phi_Folding(const int N_fold = 16)
1016 {
1017   SetOKStyle();
1018   gStyle->SetOptStat(0);
1019   gStyle->SetOptFit(1111);
1020   TVirtualFitter::SetDefaultFitter("Minuit2");
1021   gSystem->Load("libg4eval.so");
1022 
1023   TH2F * hLeakage_PHG4Particle0_e = (TH2F *) gDirectory->GetObjectChecked(
1024       "hLeakage_PHG4Particle0_e", "TH2F");
1025   TH2F * hTotal_CEMC_E_PHG4Particle0_e = (TH2F *) gDirectory->GetObjectChecked(
1026       "hTotal_CEMC_E_PHG4Particle0_e", "TH2F");
1027 
1028   TH2F * hLeakage_PHG4Particle0_e_Fold = new TH2F(
1029       "hLeakage_PHG4Particle0_e_Fold", ";Phi (rad);Leakage Ratio", 320 / N_fold,
1030       -3.14159 / N_fold, 3.14159 / N_fold, 100, 0, 0.25);
1031   TH2F * hTotal_CEMC_E_PHG4Particle0_e_Fold = new TH2F(
1032       "hTotal_CEMC_E_PHG4Particle0_e_Fold", ";Phi (rad);EMCal Energy Ratio",
1033       320 / N_fold, -3.14159 / N_fold, 3.14159 / N_fold, 100, 0, 1.2);
1034 
1035   for (int biny = 1; biny <= hLeakage_PHG4Particle0_e->GetNbinsY(); biny++)
1036     for (int bin = 1; bin <= hLeakage_PHG4Particle0_e->GetNbinsX(); bin++)
1037       {
1038         const int folded_bin = (bin - 1) % (320 / N_fold) + 1;
1039 
1040         hLeakage_PHG4Particle0_e_Fold->SetBinContent(folded_bin, biny,
1041             hLeakage_PHG4Particle0_e_Fold->GetBinContent(folded_bin, biny)
1042                 + hLeakage_PHG4Particle0_e->GetBinContent(bin, biny));
1043         hTotal_CEMC_E_PHG4Particle0_e_Fold->SetBinContent(folded_bin, biny,
1044             hTotal_CEMC_E_PHG4Particle0_e_Fold->GetBinContent(folded_bin, biny)
1045                 + hTotal_CEMC_E_PHG4Particle0_e->GetBinContent(bin, biny));
1046 
1047       }
1048 
1049   TCanvas *c1 = new TCanvas("DrawLeakage_Phi_Folding" + cuts,
1050       "DrawLeakage_Phi_Folding" + cuts, 1000, 900);
1051   c1->Divide(1, 2);
1052   int idx = 1;
1053   TPad * p;
1054 
1055   p = (TPad *) c1->cd(idx++);
1056   c1->Update();
1057   hLeakage_PHG4Particle0_e_Fold->Draw("colz");
1058 
1059   p = (TPad *) c1->cd(idx++);
1060   c1->Update();
1061   hTotal_CEMC_E_PHG4Particle0_e_Fold->Draw("colz");
1062 
1063   SaveCanvas(c1, TString("_DrawJet_") + TString(c1->GetName()), kFALSE);
1064 }
1065 
1066 void
1067 Sampling(TString infile)
1068 {
1069 
1070   TCanvas *c1 = new TCanvas("Sampling" + cuts, "Sampling" + cuts, 900, 900);
1071   c1->Divide(1, 1);
1072   int idx = 1;
1073   TPad * p = gPad;
1074 
1075   T->Draw(
1076       "Sum$(G4HIT_CEMC.light_yield):(Sum$(G4HIT_CEMC.edep) + Sum$(G4HIT_ABSORBER_CEMC.edep))",
1077       "", "*");
1078   TGraph * gh2D =
1079       (TGraph *) gPad->GetListOfPrimitives()->FindObject("Graph")->Clone(
1080           "gh2D");
1081   assert(gh2D);
1082   gh2D->SetName("gh2D");
1083   gh2D->SetMarkerStyle(kDot);
1084 
1085   p->DrawFrame(0, 0, 10, 1,
1086       "EMC Energy Deposition;Absorber+Scintillator (GeV);Scintillator (GeV)");
1087   gh2D->Draw("p");
1088 
1089   TF1 * fsampling = new TF1("fsampling", "[0]*x", 0.4, 100);
1090   gh2D->Fit(fsampling, "MR");
1091 
1092   SaveCanvas(c1,
1093       TString(_infile->GetName()) + TString("_DrawJet_")
1094           + TString(c1->GetName()), kFALSE);
1095 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!