Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 // $Id: $                                                                                             
0002 
0003 /*!
0004  * \file Draw.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 <TLatex.h>
0017 #include <TGraphErrors.h>
0018 #include <cassert>
0019 #include "SaveCanvas.C"
0020 #include "SetOKStyle.C"
0021 using namespace std;
0022 
0023 //#include "Prototype2_DSTReader.h"
0024 
0025 TCut event_sel;
0026 TString cuts;
0027 TFile * _file0 = NULL;
0028 TTree * T = NULL;
0029 
0030 class lin_res
0031 {
0032 public:
0033   TString name;
0034   TGraphErrors * linearity;
0035   TGraphErrors * resolution;
0036   TF1 * f_res;
0037 };
0038 
0039 void
0040 DrawPrototype2ShowerCalib( //
0041     const TString infile =
0042 //        "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/ShowerCalib/UIUC21.lst_EMCalCalib.root" //
0043 //        "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/ShowerCalib/THP.lst_EMCalCalib.root"//
0044 //                "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/ShowerCalib/Tilt0.lst_EMCalCalib.root"//
0045 //                "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/ShowerCalib/UpTilt5.lst_EMCalCalib.root"//
0046 //                    "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/ShowerCalib/Rot45.lst_EMCalCalib.root"//
0047 //        "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/ShowerCalib/Rot45.lst_EMCalCalib.root"//
0048 //        "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./0Degree_1Col_LightCollectionSeanStoll_CrackScan/Prototype_e-_8_SegALL_EMCalCalib.root"//
0049 //                "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./10DegreeRot_1Col_LightCollectionSeanStoll_CrackScan/Prototype_e-_8_SegALL_EMCalCalib.root"//
0050 //        "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./10DegreeRot_1Col_LightCollectionSeanStoll/Prototype_e-_ALL_SegA_ALL_EMCalCalib.root"//
0051         //        "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./10DegreeRot_1Col_LightCollectionSeanStoll/Prototype_e-_8_SegA_ALL_EMCalCalib.root"//
0052 //        "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./0Degree_1Col_LightCollectionSeanStoll_BirkConst0.151/Prototype_e-_8_SegALL_EMCalCalib.root"//
0053 //        "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./0Degree_1Col_LightCollectionSeanStoll_BirkConst0.18/Prototype_e-_8_SegALL_EMCalCalib.root"//
0054 //            "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./0Degree_1Col_LightCollectionSeanStoll/Prototype_e-_ALL_SegA_ALL_EMCalCalib.root"//
0055 //        "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./0Degree_1Col_LightCollectionSeanStoll/Prototype_pi-_8_SegALL_EMCalCalib.root"//
0056 //            "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./0Degree_1Col_LightCollectionSeanStoll/Prototype_kaon-_8_SegALL_EMCalCalib.root"//
0057 //                "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./0Degree_1Col_LightCollectionSeanStoll/Prototype_mu-_8_SegALL_EMCalCalib.root"//
0058         //        "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./45DegreeRot_1Col_LightCollectionSeanStoll/Prototype_e-_ALL_SegA_ALL_EMCalCalib.root"//
0059                 "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2017/sim/EMCalOnly_0Tilt_FlatLightCollection/Prototype_e-_ALL_SegA_ALL_EMCalCalib.root"//
0060     )
0061 {
0062 
0063   SetOKStyle();
0064   gStyle->SetOptStat(0);
0065   gStyle->SetOptFit(1111);
0066   gStyle->SetPadGridX(0);
0067   gStyle->SetPadGridY(0);
0068   TVirtualFitter::SetDefaultFitter("Minuit2");
0069 
0070   gSystem->Load("libPrototype2.so");
0071   gSystem->Load("libProto2ShowCalib.so");
0072 
0073 //  gROOT->LoadMacro("Prototype2_DSTReader.C+");
0074 
0075   if (!_file0)
0076     {
0077       TString chian_str = infile;
0078       chian_str.ReplaceAll("ALL", "*");
0079 
0080       TChain * t = new TChain("T");
0081       const int n = t->Add(chian_str);
0082 
0083       cout << "Loaded " << n << " root files with " << chian_str << endl;
0084       assert(n > 0);
0085 
0086       T = t;
0087 
0088       _file0 = new TFile;
0089       _file0->SetName(infile);
0090     }
0091 
0092   assert(_file0);
0093 
0094   event_sel = "1";
0095   cuts = "_all_data";
0096 //    event_sel = "abs(truth_z )<0.25 && abs(truth_y )<0.25";
0097 //    cuts = "_0DegreeRot_h1_v1";
0098 //  event_sel = "abs(truth_z )<1.25 && abs(truth_y )<1.25";
0099 //  cuts = "_0DegreeRot_h5_v5";
0100 //  event_sel = "abs(truth_z )<1.25 && abs(truth_y )<1.25 && abs(clus_5x5_prod.average_col -2)<=1 && abs(clus_5x5_prod.average_row -2)<=1";
0101 //  cuts = "_0DegreeRot_h5_v5_col2_row2";
0102 //  event_sel = "abs(truth_z )<.25 && abs(truth_y )<.25 && abs(clus_5x5_prod.average_col -2)<=1 && abs(clus_5x5_prod.average_row -2)<=1";
0103 //  cuts = "_0DegreeRot_h1_v1_col2_row2";
0104 //    event_sel = "abs(truth_z + 15.9)<0.25 && abs(truth_y - .1)<0.25";
0105 //    cuts = "_10DegreeRot_h1_v1";
0106 //    event_sel = "abs(truth_z + 15.9)<0.75 && abs(truth_y - .1)<0.75";
0107 //    cuts = "_10DegreeRot_h3_v3";
0108 //    event_sel = "abs(truth_z + 37)<0.25 && abs(truth_y )<0.25";
0109 //    cuts = "_45DegreeRot_h1_v1";
0110 //  event_sel = "good_data";
0111 //  cuts = "_good_data";
0112 //  event_sel = "info.beam_mom == -8 && good_data";
0113 //  cuts = "_8GeV_good_data";
0114 //  event_sel = "info.beam_mom == -8";
0115 //  cuts = "_Neg8GeV";
0116 //    event_sel = "good_data && info.hodo_h==4 && info.hodo_v==3";
0117 //    cuts = "_good_data_h4_v3";
0118   //    event_sel = "good_data && info.hodo_h>=3 && info.hodo_h<=4 && info.hodo_v>=2 && info.hodo_v<=4";
0119   //    cuts = "_good_data_h34_v234";
0120 //  event_sel = "good_data && info.hodo_h>=2 && info.hodo_h<=4 && info.hodo_v>=4 && info.hodo_v<=6";
0121 //  cuts = "_good_data_h234_v456";
0122 //  event_sel = "good_data && info.hodo_h>=3 && info.hodo_h<=3 && info.hodo_v>=5 && info.hodo_v<=5";
0123 //  cuts = "_good_data_h3_v5";
0124 //  event_sel =
0125 //      "info.beam_mom == -8 && good_temp && valid_hodo_v && valid_hodo_h&& trigger_veto_pass && info.hodo_h>=1 && info.hodo_h<=5 && info.hodo_v>=3 && info.hodo_v<=7";
0126 //  cuts = "_Neg8GeV_quality_h12345_v34567";
0127 //    event_sel =
0128 //        "info.beam_mom == -4 && good_temp && valid_hodo_v && valid_hodo_h&& trigger_veto_pass && info.hodo_h>=1 && info.hodo_h<=5 && info.hodo_v>=3 && info.hodo_v<=7 && abs(clus_5x5_prod.average_col -2)<=1 && abs(clus_5x5_prod.average_row -2)<=1";
0129 //    cuts = "_Neg4GeV_quality_h12345_v34567_col2_row2";
0130 //    event_sel =
0131 //        "info.beam_mom == -8 && good_temp && valid_hodo_v && valid_hodo_h&& trigger_veto_pass && info.hodo_h>=1 && info.hodo_h<=5 && info.hodo_v>=3 && info.hodo_v<=7 && abs(clus_5x5_prod.average_col -2)<=1 && abs(clus_5x5_prod.average_row -2)<=1";
0132 //    cuts = "_Neg8GeV_quality_h12345_v34567_col2_row2";
0133 //    event_sel =
0134 //        "info.beam_mom == -12 && good_temp && valid_hodo_v && valid_hodo_h&& trigger_veto_pass && info.hodo_h>=1 && info.hodo_h<=5 && info.hodo_v>=3 && info.hodo_v<=7 && abs(clus_5x5_prod.average_col -2)<=1 && abs(clus_5x5_prod.average_row -2)<=1";
0135 //    cuts = "_Neg12GeV_quality_h12345_v34567_col2_row2";
0136 //  event_sel =
0137 //      "info.beam_mom == -4 && good_temp && valid_hodo_v && valid_hodo_h&& trigger_veto_pass && info.hodo_h==3 && info.hodo_v==5 && abs(clus_5x5_prod.average_col -2)<=1 && abs(clus_5x5_prod.average_row -2)<=1";
0138 //  cuts = "_Neg4GeV_quality_h3_v5_col2_row2";
0139 //  event_sel =
0140 //      "info.beam_mom == -8 && good_temp && valid_hodo_v && valid_hodo_h&& trigger_veto_pass && info.hodo_h==3 && info.hodo_v==5 && abs(clus_5x5_prod.average_col -2)<=1 && abs(clus_5x5_prod.average_row -2)<=1";
0141 //  cuts = "_Neg8GeV_quality_h3_v5_col2_row2";
0142 //  event_sel =
0143 //      "info.beam_mom == -12 && good_temp && valid_hodo_v && valid_hodo_h&& trigger_veto_pass && info.hodo_h==3 && info.hodo_v==5 && abs(clus_5x5_prod.average_col -2)<=1 && abs(clus_5x5_prod.average_row -2)<=1";
0144 //  cuts = "_Neg12GeV_quality_h3_v5_col2_row2";
0145 //    event_sel = // UpTilt5
0146 //        "info.beam_mom == -8 && good_temp && valid_hodo_v && valid_hodo_h&& trigger_veto_pass && info.hodo_h==4 && info.hodo_v==6 && abs(clus_5x5_prod.average_col -2)<=1 && abs(clus_5x5_prod.average_row -2)<=1";
0147 //    cuts = "_Neg8GeV_quality_h4_v6_col2_row2";
0148 //    event_sel = "good_data && info.hodo_h>=4 && info.hodo_h<=6 && info.hodo_v>=2 && info.hodo_v<=4";
0149 //    cuts = "_good_data_h456_v234";
0150 //  event_sel = "good_data && info.hodo_h==5 && info.hodo_v==3";
0151 //  cuts = "_good_data_h5_v3";
0152 //  event_sel = "info.beam_mom == -8 && good_data && info.hodo_h==5 && info.hodo_v==3";
0153 //  cuts = "_Neg8GeV_good_data_h5_v3";
0154 //      event_sel = "good_data && info.hodo_h==5 && (info.hodo_v==4 || info.hodo_v==5)";
0155 //      cuts = "_good_data_h5_v54";
0156 
0157 //  event_sel = "good_data && info.hodo_h==5 && info.hodo_v==5";
0158 //  cuts = "_good_data_h5_v5";
0159 //  event_sel = "good_data &&  info.hodo_v==5";
0160 //  cuts = "_good_data_hall_v5";
0161 //  event_sel = "good_data && ( info.hodo_v==5 ||  info.hodo_v==4)";
0162 //  cuts = "_good_data_hall_v45";
0163 
0164       T->SetAlias("SimEnergyScale","1*1");
0165 //    // based on /phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/ShowerCalib/UIUC21.lst_EMCalCalib.root_DrawPrototype2ShowerCalib_LineShapeData_Neg8GeV_good_data_h5_v3.svg
0166 //    T->SetAlias("SimEnergyScale","8.74635e+00/7.60551");
0167   //  // based on /phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/ShowerCalib/Tilt0.lst_EMCalCalib.root_DrawPrototype2ShowerCalib_LineShapeData_Neg8GeV_quality_h3_v5_col2_row2.root_DrawPrototype2ShowerCalib_SumLineShapeCompare_Electron_8GeV_QGSP_BERT_HP.root
0168 //    T->SetAlias("SimEnergyScale","8.88178/8.16125e+00");
0169   // Tilt0.lst <-> QGSP_BERT_HP Birk 0.151
0170 //    T->SetAlias("SimEnergyScale","8.88178/8.01845");
0171   // Tilt0.lst <-> QGSP_BERT_HP Birk 0.18
0172 //    T->SetAlias("SimEnergyScale","8.88178/7.94838");
0173 
0174     //  Dummy
0175 //    T->SetAlias("SimEnergyScale","1*1");
0176 
0177   cout << "Build event selection of " << (const char *) event_sel << endl;
0178 
0179   T->Draw(">>EventList", event_sel);
0180   TEventList * elist = gDirectory->GetObjectChecked("EventList", "TEventList");
0181   cout << elist->GetN() << " / " << T->GetEntriesFast() << " events selected"
0182       << endl;
0183 
0184   T->SetEventList(elist);
0185 
0186 //  // data stuff
0187 //  PositionDependenceData("clus_5x5_prod.sum_E");
0188 //  PositionDependenceData("clus_5x5_recalib.sum_E");
0189 //  HodoscopeCheck();
0190 //    LineShapeData("abs(info.C2_sum)<100",  "(info.C2_sum)>600 && (info.C2_sum)<1300"); // 4 GeV
0191 //    LineShapeData("abs(info.C2_sum)<100",  "(info.C2_sum)>500 && (info.C2_sum)<1300"); // 8 GeV
0192 //  LineShapeData("abs(info.C2_sum)<100",  "(info.C2_sum)>200 && (info.C2_sum)<1300"); // 12 GeV
0193 
0194 //  Get_Res_linear_Summmary();
0195 
0196   // simulation stuff
0197 //  SimPositionCheck(-0); // 0 degree tilted
0198 //  LineShapeSim();
0199 
0200   PositionDependenceSim("clus_5x5_prod.sum_E", -0, 5); // 0 degree tilted
0201 //  SimPositionCheck(-15); // 10 degree tilted
0202 //  PositionDependenceSim("clus_5x5_prod.sum_E", -15, 5); // 10 degree tilted
0203 //    SimPositionCheck(-40+3); // 45 degree tilted
0204   Get_Res_linear_Summmary_Sim();
0205 }
0206 
0207 void
0208 PositionDependenceData(TString sTOWER = "clus_5x5_prod.sum_E",
0209     const double z_shift = 0, const int n_div = 1)
0210 {
0211   TH3F * EnergySum_LG3 =
0212       new TH3F("EnergySum_LG3",
0213           ";Horizontal Hodoscope (5 mm);Vertical Hodoscope (5 mm);5x5 Cluster Energy (GeV)", //
0214           8, -.5, 7.5, //
0215           8, -.5, 7.5, //
0216           200, 1, 20);
0217 
0218   T->Draw(sTOWER + ":7-hodo_v:hodo_h>>EnergySum_LG3", "", "goff");
0219 
0220   TProfile2D * EnergySum_LG3_xy = EnergySum_LG3->Project3DProfile("yx");
0221   TH2 * EnergySum_LG3_zx = EnergySum_LG3->Project3D("zx");
0222   TH2 * EnergySum_LG3_zy = EnergySum_LG3->Project3D("zy");
0223 
0224   TGraphErrors * ge_EnergySum_LG3_zx = FitProfile(EnergySum_LG3_zx);
0225   TGraphErrors * ge_EnergySum_LG3_zy = FitProfile(EnergySum_LG3_zy);
0226 
0227   TText * t;
0228   TCanvas *c1 = new TCanvas(
0229       TString(Form("EMCDistributionVSBeam_SUM_NDiv%d_", n_div)) + sTOWER + cuts,
0230       TString(Form("EMCDistributionVSBeam_SUM_NDiv%d_", n_div)) + sTOWER + cuts,
0231       1000, 960);
0232   c1->Divide(2, 2);
0233   int idx = 1;
0234   TPad * p;
0235 
0236   p = (TPad *) c1->cd(idx++);
0237   c1->Update();
0238 //  p->SetLogy();
0239   p->SetGridx(0);
0240   p->SetGridy(0);
0241 
0242   EnergySum_LG3_xy->Draw("colz");
0243   EnergySum_LG3_xy->SetTitle(
0244       "Position scan;Horizontal Hodoscope (5 mm);Vertical Hodoscope (5 mm)");
0245 
0246   p = (TPad *) c1->cd(idx++);
0247   c1->Update();
0248 //  p->SetLogy();
0249   p->SetGridx(0);
0250   p->SetGridy(0);
0251 
0252   EnergySum_LG3_zx->Draw("colz");
0253   EnergySum_LG3_zx->SetTitle(
0254       "Position scan;Horizontal Hodoscope (5 mm);5x5 Cluster Energy (GeV)");
0255 
0256   ge_EnergySum_LG3_zx->SetLineWidth(2);
0257   ge_EnergySum_LG3_zx->SetMarkerStyle(kFullCircle);
0258   ge_EnergySum_LG3_zx->Draw("pe");
0259 
0260   p = (TPad *) c1->cd(idx++);
0261   c1->Update();
0262 //  p->SetLogy();
0263   p->SetGridx(0);
0264   p->SetGridy(0);
0265 
0266   EnergySum_LG3_zy->Draw("colz");
0267   EnergySum_LG3_zy->SetTitle(
0268       "Position scan;Vertical Hodoscope (5 mm);5x5 Cluster Energy (GeV)");
0269 
0270   ge_EnergySum_LG3_zy->SetLineWidth(2);
0271   ge_EnergySum_LG3_zy->SetMarkerStyle(kFullCircle);
0272   ge_EnergySum_LG3_zy->Draw("pe");
0273 
0274   p = (TPad *) c1->cd(idx++);
0275   c1->Update();
0276 //  p->SetLogy();
0277   p->SetGridx(0);
0278   p->SetGridy(0);
0279 
0280   TH1 * h = (TH1 *) EnergySum_LG3->ProjectionZ();
0281 
0282   TF1 * fgaus = new TF1("fgaus_LG", "gaus", 0, 100);
0283   fgaus->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
0284       h->GetMean() + 2 * h->GetRMS());
0285   h->Fit(fgaus, "M");
0286 
0287   h->Sumw2();
0288   h->GetXaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
0289       h->GetMean() + 4 * h->GetRMS());
0290   EnergySum_LG3_zx->GetYaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
0291       h->GetMean() + 4 * h->GetRMS());
0292   EnergySum_LG3_zy->GetYaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
0293       h->GetMean() + 4 * h->GetRMS());
0294 
0295   h->SetLineWidth(2);
0296   h->SetMarkerStyle(kFullCircle);
0297 
0298   h->SetTitle(
0299       Form("#DeltaE/<E> = %.1f%%",
0300           100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));
0301 
0302   SaveCanvas(c1,
0303       TString(_file0->GetName()) + TString("_DrawPrototype2ShowerCalib_")
0304           + TString(c1->GetName()), kTRUE);
0305 }
0306 
0307 void
0308 PositionDependenceSim(TString sTOWER = "clus_5x5_prod.sum_E",
0309     const double z_shift = 0, const int n_div = 1)
0310 {
0311   TH3F * EnergySum_LG3 =
0312       new TH3F("EnergySum_LG3",
0313           ";Beam Horizontal Pos (cm);Beam Vertical Pos (cm);5x5 Cluster Energy (GeV)", //
0314           20 * n_div, z_shift - 5, z_shift + 5, //
0315           20 * n_div, -5, 5, //
0316           200, 1, 20);
0317 
0318   T->Draw(sTOWER + ":info.truth_y:info.truth_z>>EnergySum_LG3", "", "goff");
0319 
0320   TProfile2D * EnergySum_LG3_xy = EnergySum_LG3->Project3DProfile("yx");
0321   TH2 * EnergySum_LG3_zx = EnergySum_LG3->Project3D("zx");
0322   TH2 * EnergySum_LG3_zy = EnergySum_LG3->Project3D("zy");
0323 
0324   TGraphErrors * ge_EnergySum_LG3_zx = FitProfile(EnergySum_LG3_zx);
0325   TGraphErrors * ge_EnergySum_LG3_zy = FitProfile(EnergySum_LG3_zy);
0326 
0327   TText * t;
0328   TCanvas *c1 = new TCanvas(
0329       TString(Form("EMCDistributionVSBeam_SUM_NDiv%d_", n_div)) + sTOWER + cuts,
0330       TString(Form("EMCDistributionVSBeam_SUM_NDiv%d_", n_div)) + sTOWER + cuts,
0331       1000, 960);
0332   c1->Divide(2, 2);
0333   int idx = 1;
0334   TPad * p;
0335 
0336   p = (TPad *) c1->cd(idx++);
0337   c1->Update();
0338 //  p->SetLogy();
0339   p->SetGridx(0);
0340   p->SetGridy(0);
0341 
0342   EnergySum_LG3_xy->Draw("colz");
0343   EnergySum_LG3_xy->SetTitle(
0344       "Position scan;Beam Horizontal Pos (cm);Beam Vertical Pos (cm)");
0345 
0346   p = (TPad *) c1->cd(idx++);
0347   c1->Update();
0348 //  p->SetLogy();
0349   p->SetGridx(0);
0350   p->SetGridy(0);
0351 
0352   EnergySum_LG3_zx->Draw("colz");
0353   EnergySum_LG3_zx->SetTitle(
0354       "Position scan;Beam Horizontal Pos (cm);5x5 Cluster Energy (GeV)");
0355 
0356   ge_EnergySum_LG3_zx->SetLineWidth(2);
0357   ge_EnergySum_LG3_zx->SetMarkerStyle(kFullCircle);
0358   ge_EnergySum_LG3_zx->Draw("pe");
0359 
0360   p = (TPad *) c1->cd(idx++);
0361   c1->Update();
0362 //  p->SetLogy();
0363   p->SetGridx(0);
0364   p->SetGridy(0);
0365 
0366   EnergySum_LG3_zy->Draw("colz");
0367   EnergySum_LG3_zy->SetTitle(
0368       "Position scan;Beam Vertical Pos (cm);5x5 Cluster Energy (GeV)");
0369 
0370   ge_EnergySum_LG3_zy->SetLineWidth(2);
0371   ge_EnergySum_LG3_zy->SetMarkerStyle(kFullCircle);
0372   ge_EnergySum_LG3_zy->Draw("pe");
0373 
0374   p = (TPad *) c1->cd(idx++);
0375   c1->Update();
0376 //  p->SetLogy();
0377   p->SetGridx(0);
0378   p->SetGridy(0);
0379 
0380   TH1 * h = (TH1 *) EnergySum_LG3->ProjectionZ();
0381 
0382   TF1 * fgaus = new TF1("fgaus_LG", "gaus", 0, 100);
0383   fgaus->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
0384       h->GetMean() + 2 * h->GetRMS());
0385   h->Fit(fgaus, "M");
0386 
0387   h->Sumw2();
0388   h->GetXaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
0389       h->GetMean() + 4 * h->GetRMS());
0390   EnergySum_LG3_zx->GetYaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
0391       h->GetMean() + 4 * h->GetRMS());
0392   EnergySum_LG3_zy->GetYaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
0393       h->GetMean() + 4 * h->GetRMS());
0394 
0395   h->SetLineWidth(2);
0396   h->SetMarkerStyle(kFullCircle);
0397 
0398   h->SetTitle(
0399       Form("#DeltaE/<E> = %.1f%%",
0400           100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));
0401 
0402   SaveCanvas(c1,
0403       TString(_file0->GetName()) + TString("_DrawPrototype2ShowerCalib_")
0404           + TString(c1->GetName()), kTRUE);
0405 }
0406 
0407 void
0408 LineShapeData(TCut c2_h = "abs(info.C2_sum)<100", TCut c2_e =
0409     "(info.C2_sum)>500 && (info.C2_sum)<1300")
0410 {
0411 
0412   vector<double> mom;
0413 
0414   TText * t;
0415   TCanvas *c1 = new TCanvas("LineShapeData" + cuts, "LineShapeData" + cuts,
0416       1800, 950);
0417 
0418   c1->Divide(2, 2);
0419   int idx = 1;
0420   TPad * p;
0421 
0422   p = (TPad *) c1->cd(idx++);
0423   c1->Update();
0424 //  p->SetGridx(1);
0425 //  p->SetGridy(1);
0426   p->SetLogy();
0427 
0428   T->Draw("info.C2_sum>>h_c2_sum(300,-500,2500)");
0429   h_c2_sum->SetTitle("Cherenkov Checks;Sum C2 (ADC)");
0430   T->Draw("info.C2_sum>>h_c2_h(300,-500,2500)", c2_h, "same");
0431   T->Draw("info.C2_sum>>h_c2_e(300,-500,2500)", c2_e, "same");
0432 
0433   h_c2_h->SetLineColor(kBlue);
0434   h_c2_e->SetLineColor(kRed);
0435 
0436   h_c2_sum->SetLineWidth(2);
0437   h_c2_h->SetLineWidth(2);
0438   h_c2_e->SetLineWidth(2);
0439 
0440   p = (TPad *) c1->cd(idx++);
0441   c1->Update();
0442   p->SetLogy();
0443 
0444   T->Draw("clus_5x5_recalib.sum_E>>h_5x5sum_c2_sum(170,-1,16)");
0445   h_5x5sum_c2_sum->SetTitle(
0446       "Cluster spectrum decomposition;5x5 cluster energy (GeV)");
0447   T->Draw("clus_5x5_recalib.sum_E>>h_5x5sum_c2_h(170,-1,16)", c2_h, "same");
0448   T->Draw("clus_5x5_recalib.sum_E>>h_5x5sum_c2_rej_h(170,-1,16)", !c2_h,
0449       "same");
0450   T->Draw("clus_5x5_recalib.sum_E>>h_5x5sum_c2_e(170,-1,16)", c2_e, "same");
0451 
0452   h_5x5sum_c2_h->SetLineColor(kBlue);
0453   h_5x5sum_c2_rej_h->SetLineColor(kMagenta);
0454   h_5x5sum_c2_e->SetLineColor(kRed);
0455 
0456   h_5x5sum_c2_sum->SetLineWidth(2);
0457   h_5x5sum_c2_h->SetLineWidth(2);
0458   h_5x5sum_c2_rej_h->SetLineWidth(2);
0459   h_5x5sum_c2_e->SetLineWidth(2);
0460 
0461   p = (TPad *) c1->cd(idx++);
0462   c1->Update();
0463   p->SetLogy();
0464 
0465   TH1 * h_5x5sum_c2_h2 = h_5x5sum_c2_h->DrawClone();
0466   h_5x5sum_c2_h2->Sumw2();
0467   h_5x5sum_c2_h2->SetMarkerColor(kBlue);
0468   h_5x5sum_c2_h2->SetMarkerStyle(kFullCircle);
0469   h_5x5sum_c2_h2->SetTitle(";5x5 cluster energy (GeV)");
0470 
0471   TH1 * h_5x5sum_c2_e2 = h_5x5sum_c2_e->DrawClone("same");
0472   h_5x5sum_c2_e2->Sumw2();
0473   h_5x5sum_c2_e2->SetMarkerColor(kRed);
0474   h_5x5sum_c2_e2->SetMarkerStyle(kFullCircle);
0475   h_5x5sum_c2_e2->SetTitle(";5x5 cluster energy (GeV)");
0476 
0477   p = (TPad *) c1->cd(idx++);
0478   c1->Update();
0479   p->SetLogy();
0480 
0481   TH1F * h_5x5sum_c2_h3 = h_5x5sum_c2_h->DrawClone();
0482   h_5x5sum_c2_h3->SetName("h_5x5sum_c2_h3");
0483   h_5x5sum_c2_h3->Sumw2();
0484   h_5x5sum_c2_h3->Scale(1. / h_5x5sum_c2_h3->GetSum());
0485   h_5x5sum_c2_h3->SetMarkerColor(kBlue);
0486   h_5x5sum_c2_h3->SetMarkerStyle(kFullCircle);
0487   h_5x5sum_c2_h3->SetTitle(";5x5 cluster energy (GeV);Probability / bin");
0488 
0489   SaveCanvas(c1,
0490       TString(_file0->GetName()) + "_DrawPrototype2ShowerCalib_"
0491           + TString(c1->GetName()), kTRUE);
0492 }
0493 
0494 void
0495 LineShapeSim()
0496 {
0497   vector<double> mom;
0498 
0499   TText * t;
0500   TCanvas *c1 = new TCanvas("LineShapeSim" + cuts, "LineShapeSim" + cuts, 1000,
0501       650);
0502 
0503 //  c1->Divide(2,2);
0504   int idx = 1;
0505   TPad * p;
0506 
0507   p = (TPad *) c1->cd(idx++);
0508   c1->Update();
0509   p->SetLogy();
0510 
0511   T->Draw("clus_5x5_prod.sum_E*SimEnergyScale>>h_5x5sum_c2_sum(170,-1,16)");
0512   h_5x5sum_c2_sum->SetTitle(
0513       "Cluster spectrum decomposition;5x5 cluster energy (GeV)");
0514   h_5x5sum_c2_sum->SetLineWidth(2);
0515 
0516   SaveCanvas(c1,
0517       TString(_file0->GetName()) + "_DrawPrototype2ShowerCalib_"
0518           + TString(c1->GetName()), kTRUE);
0519 
0520 }
0521 
0522 void
0523 HodoscopeCheck()
0524 {
0525   vector<double> mom;
0526 
0527   TText * t;
0528   TCanvas *c1 = new TCanvas("HodoscopeCheck" + cuts, "HodoscopeCheck" + cuts,
0529       1300, 950);
0530 
0531   c1->Divide(2, 1);
0532   int idx = 1;
0533   TPad * p;
0534 
0535   p = (TPad *) c1->cd(idx++);
0536   c1->Update();
0537   p->SetGridx(1);
0538   p->SetGridy(1);
0539   p->SetLogz();
0540 
0541   T->Draw("clus_5x5_prod.average_col:hodo_h>>h2_h(8,-.5,7.5,160,-.5,7.5)",
0542       "good_data", "colz");
0543   h2_h->SetTitle(
0544       "Horizontal hodoscope check;Horizontal Hodoscope;5x5 cluster mean col");
0545 
0546   p = (TPad *) c1->cd(idx++);
0547   c1->Update();
0548   p->SetGridx(1);
0549   p->SetGridy(1);
0550   p->SetLogz();
0551 
0552   T->Draw("clus_5x5_prod.average_row:hodo_v>>h2_v(8,-.5,7.5,160,-.5,7.5)",
0553       "good_data", "colz");
0554   h2_v->SetTitle(
0555       "Vertical hodoscope check;Vertical Hodoscope;5x5 cluster mean row");
0556 
0557   SaveCanvas(c1,
0558       TString(_file0->GetName()) + "_DrawPrototype2ShowerCalib_"
0559           + TString(c1->GetName()), kTRUE);
0560 
0561   return mom;
0562 }
0563 
0564 void
0565 SimPositionCheck(const double shift_z = 0)
0566 {
0567   vector<double> mom;
0568 
0569   TText * t;
0570   TCanvas *c1 = new TCanvas("SimPositionCheck" + cuts,
0571       "SimPositionCheck" + cuts, 1300, 950);
0572 
0573   c1->Divide(2, 1);
0574   int idx = 1;
0575   TPad * p;
0576 
0577   p = (TPad *) c1->cd(idx++);
0578   c1->Update();
0579   p->SetGridx(1);
0580   p->SetGridy(1);
0581   p->SetLogz();
0582 
0583   T->Draw(
0584       Form("clus_5x5_prod.average_col:truth_z>>h2_h(30,%f,%f,160,-.5,7.5)",
0585           shift_z - 1.5, shift_z + 1.5), "1", "colz");
0586   h2_h->SetTitle(
0587       "Horizontal hodoscope check;Horizontal beam pos;5x5 cluster mean col");
0588 
0589   p = (TPad *) c1->cd(idx++);
0590   c1->Update();
0591   p->SetGridx(1);
0592   p->SetGridy(1);
0593   p->SetLogz();
0594 
0595   T->Draw("clus_5x5_prod.average_row:truth_y>>h2_v(30,-1.5,1.5,160,-.5,7.5)",
0596       "1", "colz");
0597   h2_v->SetTitle(
0598       "Vertical hodoscope check;Vertical beam pos;5x5 cluster mean row");
0599 
0600   SaveCanvas(c1,
0601       TString(_file0->GetName()) + "_DrawPrototype2ShowerCalib_"
0602           + TString(c1->GetName()), kTRUE);
0603 
0604   return;
0605 }
0606 
0607 void
0608 Get_Res_linear_Summmary()
0609 {
0610 
0611   vector<double> beam_mom(GetBeamMom());
0612 
0613 //  return;
0614 
0615   lin_res ges_clus_5x5_prod = GetResolution("clus_5x5_prod", beam_mom, kBlue);
0616   lin_res ges_clus_3x3_prod = GetResolution("clus_3x3_prod", beam_mom,
0617       kBlue + 3);
0618   lin_res ges_clus_5x5_temp = GetResolution("clus_5x5_temp", beam_mom,
0619       kRed - 2);
0620   lin_res ges_clus_5x5_recalib = GetResolution("clus_5x5_recalib", beam_mom,
0621       kRed + 3);
0622 
0623   TCanvas *c1 = new TCanvas(Form("Res_linear") + cuts,
0624       Form("Res_linear") + cuts, 1300, 600);
0625   c1->Divide(2, 1);
0626   int idx = 1;
0627   TPad * p;
0628 
0629   p = (TPad *) c1->cd(idx++);
0630   c1->Update();
0631 //  p->SetLogy();
0632 
0633   TLegend* leg = new TLegend(.15, .7, .6, .85);
0634 
0635   p->DrawFrame(0, 0, 25, 25,
0636       Form("Electron Linearity;Input energy (GeV);Measured Energy (GeV)"));
0637   TLine * l = new TLine(0, 0, 25, 25);
0638   l->SetLineColor(kGray);
0639 //  l->Draw();
0640 
0641   TF1 * f_calo_l_sim = new TF1("f_calo_l", "pol2", 0.5, 25);
0642 //  f_calo_l_sim->SetParameters(-0.03389, 0.9666, -0.0002822);
0643   f_calo_l_sim->SetParameters(-0., 1, -0.);
0644 //  f_calo_l_sim->SetLineWidth(1);
0645   f_calo_l_sim->SetLineColor(kGreen + 2);
0646   f_calo_l_sim->SetLineWidth(3);
0647 
0648   f_calo_l_sim->Draw("same");
0649   ges_clus_5x5_prod.linearity->Draw("p");
0650   ges_clus_3x3_prod.linearity->Draw("p");
0651   ges_clus_5x5_temp.linearity->Draw("p");
0652   ges_clus_5x5_recalib.linearity->Draw("p");
0653 //  ge_linear->Fit(f_calo_l, "RM0");
0654 //  f_calo_l->Draw("same");
0655 
0656   leg->AddEntry(ges_clus_5x5_prod.linearity, ges_clus_5x5_prod.name, "ep");
0657   leg->AddEntry(ges_clus_3x3_prod.linearity, ges_clus_3x3_prod.name, "ep");
0658   leg->AddEntry(ges_clus_5x5_temp.linearity, ges_clus_5x5_temp.name, "ep");
0659   leg->AddEntry(ges_clus_5x5_recalib.linearity, "clus_5x5_recalib", "ep");
0660   leg->AddEntry(f_calo_l_sim, "Unity", "l");
0661   leg->Draw();
0662 
0663   p = (TPad *) c1->cd(idx++);
0664   c1->Update();
0665 //  p->SetLogy();
0666   p->SetGridx(0);
0667   p->SetGridy(0);
0668 
0669   TF1 * f_calo_sim = new TF1("f_calo_sim", "sqrt([0]*[0]+[1]*[1]/x)/100", 0.5,
0670       25);
0671   f_calo_sim->SetParameters(2.4, 11.8);
0672   f_calo_sim->SetLineWidth(3);
0673   f_calo_sim->SetLineColor(kGreen + 2);
0674 
0675   TH1 * hframe = p->DrawFrame(0, 0, 25, 0.3,
0676       Form("Resolution;Input energy (GeV);#DeltaE/<E>"));
0677 
0678   TLegend* leg = new TLegend(.2, .6, .85, .9);
0679 
0680   ges_clus_5x5_prod.f_res->Draw("same");
0681   ges_clus_5x5_prod.resolution->Draw("ep");
0682   ges_clus_3x3_prod.f_res->Draw("same");
0683   ges_clus_3x3_prod.resolution->Draw("ep");
0684   ges_clus_5x5_temp.f_res->Draw("same");
0685   ges_clus_5x5_temp.resolution->Draw("ep");
0686   ges_clus_5x5_recalib.f_res->Draw("same");
0687   ges_clus_5x5_recalib.resolution->Draw("ep");
0688   f_calo_sim->Draw("same");
0689 
0690   leg->AddEntry(ges_clus_5x5_prod.resolution, ges_clus_5x5_prod.name, "ep");
0691   leg->AddEntry(ges_clus_5x5_prod.f_res,
0692       Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
0693           ges_clus_5x5_prod.f_res->GetParameter(0),
0694           ges_clus_5x5_prod.f_res->GetParameter(1)), "l");
0695 
0696   leg->AddEntry(ges_clus_3x3_prod.resolution, ges_clus_3x3_prod.name, "ep");
0697   leg->AddEntry(ges_clus_3x3_prod.f_res,
0698       Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
0699           ges_clus_3x3_prod.f_res->GetParameter(0),
0700           ges_clus_3x3_prod.f_res->GetParameter(1)), "l");
0701 
0702   leg->AddEntry(ges_clus_5x5_temp.resolution, ges_clus_5x5_temp.name, "ep");
0703   leg->AddEntry(ges_clus_5x5_temp.f_res,
0704       Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
0705           ges_clus_5x5_temp.f_res->GetParameter(0),
0706           ges_clus_5x5_temp.f_res->GetParameter(1)), "l");
0707 
0708   leg->AddEntry(ges_clus_5x5_recalib.resolution, "clus_5x5_recalib", "ep");
0709   leg->AddEntry(ges_clus_5x5_recalib.f_res,
0710       Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
0711           ges_clus_5x5_recalib.f_res->GetParameter(0),
0712           ges_clus_5x5_recalib.f_res->GetParameter(1)), "l");
0713 //  leg->AddEntry(new TH1(), "", "l");
0714 //  leg->AddEntry((TObject*) 0, " ", "");
0715 
0716   leg->Draw();
0717 
0718   TLegend* leg = new TLegend(.2, .1, .85, .3);
0719 
0720   leg->AddEntry(f_calo_sim,
0721       Form("Prelim. Sim., #DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
0722           f_calo_sim->GetParameter(0), f_calo_sim->GetParameter(1)), "l");
0723   leg->Draw();
0724 
0725   hframe->SetTitle("Electron Resolution");
0726 
0727   SaveCanvas(c1,
0728       TString(_file0->GetName()) + "_DrawPrototype2ShowerCalib_"
0729           + TString(c1->GetName()), kTRUE);
0730 }
0731 
0732 void
0733 Get_Res_linear_Summmary_Sim()
0734 {
0735 
0736   vector<double> beam_mom(GetBeamMom());
0737 
0738 //  return;
0739 
0740   lin_res ges_clus_5x5_prod = GetResolution("clus_5x5_prod", beam_mom, kBlue);
0741   lin_res ges_clus_3x3_prod = GetResolution("clus_3x3_prod", beam_mom,
0742       kBlue + 3);
0743 
0744   TCanvas *c1 = new TCanvas(Form("Res_linear") + cuts,
0745       Form("Res_linear") + cuts, 1300, 600);
0746   c1->Divide(2, 1);
0747   int idx = 1;
0748   TPad * p;
0749 
0750   p = (TPad *) c1->cd(idx++);
0751   c1->Update();
0752 //  p->SetLogy();
0753 
0754   TLegend* leg = new TLegend(.15, .7, .6, .85);
0755 
0756   p->DrawFrame(0, 0, 25, 25,
0757       Form("Electron Linearity;Input energy (GeV);Measured Energy (GeV)"));
0758   TLine * l = new TLine(0, 0, 25, 25);
0759   l->SetLineColor(kGray);
0760 //  l->Draw();
0761 
0762   TF1 * f_calo_l_sim = new TF1("f_calo_l", "pol2", 0.5, 25);
0763 //  f_calo_l_sim->SetParameters(-0.03389, 0.9666, -0.0002822);
0764   f_calo_l_sim->SetParameters(-0., 1, -0.);
0765 //  f_calo_l_sim->SetLineWidth(1);
0766   f_calo_l_sim->SetLineColor(kGreen + 2);
0767   f_calo_l_sim->SetLineWidth(3);
0768 
0769   f_calo_l_sim->Draw("same");
0770   ges_clus_5x5_prod.linearity->Draw("p");
0771   ges_clus_3x3_prod.linearity->Draw("p");
0772 //  ge_linear->Fit(f_calo_l, "RM0");
0773 //  f_calo_l->Draw("same");
0774 
0775   leg->AddEntry(ges_clus_5x5_prod.linearity, ges_clus_5x5_prod.name, "ep");
0776   leg->AddEntry(ges_clus_3x3_prod.linearity, ges_clus_3x3_prod.name, "ep");
0777   leg->AddEntry(f_calo_l_sim, "Unity", "l");
0778   leg->Draw();
0779 
0780   p = (TPad *) c1->cd(idx++);
0781   c1->Update();
0782 //  p->SetLogy();
0783   p->SetGridx(0);
0784   p->SetGridy(0);
0785 
0786   TF1 * f_calo_sim = new TF1("f_calo_sim", "sqrt([0]*[0]+[1]*[1]/x)/100", 0.5,
0787       25);
0788   f_calo_sim->SetParameters(2.4, 11.8);
0789   f_calo_sim->SetLineWidth(3);
0790   f_calo_sim->SetLineColor(kGreen + 2);
0791 
0792   TH1 * hframe = p->DrawFrame(0, 0, 25, 0.3,
0793       Form("Resolution;Input energy (GeV);#DeltaE/<E>"));
0794 
0795   TLegend* leg = new TLegend(.2, .6, .85, .9);
0796 
0797   ges_clus_5x5_prod.f_res->Draw("same");
0798   ges_clus_5x5_prod.resolution->Draw("ep");
0799   ges_clus_3x3_prod.f_res->Draw("same");
0800   ges_clus_3x3_prod.resolution->Draw("ep");
0801   f_calo_sim->Draw("same");
0802 
0803   leg->AddEntry(ges_clus_5x5_prod.resolution, ges_clus_5x5_prod.name, "ep");
0804   leg->AddEntry(ges_clus_5x5_prod.f_res,
0805       Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
0806           ges_clus_5x5_prod.f_res->GetParameter(0),
0807           ges_clus_5x5_prod.f_res->GetParameter(1)), "l");
0808 
0809   leg->AddEntry(ges_clus_3x3_prod.resolution, ges_clus_3x3_prod.name, "ep");
0810   leg->AddEntry(ges_clus_3x3_prod.f_res,
0811       Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
0812           ges_clus_3x3_prod.f_res->GetParameter(0),
0813           ges_clus_3x3_prod.f_res->GetParameter(1)), "l");
0814 
0815 //  leg->AddEntry(new TH1(), "", "l");
0816 //  leg->AddEntry((TObject*) 0, " ", "");
0817 
0818   leg->Draw();
0819 
0820   TLegend* leg = new TLegend(.2, .1, .85, .3);
0821 
0822   leg->AddEntry(f_calo_sim,
0823       Form("Prelim. Sim., #DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
0824           f_calo_sim->GetParameter(0), f_calo_sim->GetParameter(1)), "l");
0825   leg->Draw();
0826 
0827   hframe->SetTitle("Electron Resolution");
0828 
0829   SaveCanvas(c1,
0830       TString(_file0->GetName()) + "_DrawPrototype2ShowerCalib_"
0831           + TString(c1->GetName()), kTRUE);
0832 }
0833 
0834 vector<double>
0835 GetBeamMom()
0836 {
0837   vector<double> mom;
0838 
0839   TH1F * hbeam_mom = new TH1F("hbeam_mom", ";beam momentum (GeV)", 32, .5,
0840       32.5);
0841 
0842   TText * t;
0843   TCanvas *c1 = new TCanvas("GetBeamMom" + cuts, "GetBeamMom" + cuts, 1800,
0844       900);
0845 
0846   T->Draw("abs(info.beam_mom)>>hbeam_mom");
0847 
0848   for (int bin = 1; bin < hbeam_mom->GetNbinsX(); bin++)
0849     {
0850       if (hbeam_mom->GetBinContent(bin) > 50)
0851         {
0852           const double momentum = hbeam_mom->GetBinCenter(bin);
0853 
0854           if (momentum == 1 || momentum == 2 || momentum == 3 || momentum == 4
0855               || momentum == 6 || momentum == 8 || momentum == 12
0856               || momentum == 16 || momentum == 24 || momentum == 32)
0857             {
0858               mom.push_back(momentum);
0859 
0860               cout << "GetBeamMom - " << momentum << " GeV for "
0861                   << hbeam_mom->GetBinContent(bin) << " event" << endl;
0862             }
0863         }
0864     }
0865 
0866   SaveCanvas(c1,
0867       TString(_file0->GetName()) + "_DrawPrototype2ShowerCalib_"
0868           + TString(c1->GetName()), kTRUE);
0869 
0870   return mom;
0871 }
0872 
0873 lin_res
0874 GetResolution(TString cluster_name, vector<double> beam_mom, Color_t col)
0875 {
0876   vector<double> mean;
0877   vector<double> mean_err;
0878   vector<double> res;
0879   vector<double> res_err;
0880 
0881   TCanvas *c1 = new TCanvas("GetResolution_LineShape_" + cluster_name + cuts,
0882       "GetResolution_LineShape_" + cluster_name + cuts, 1800, 900);
0883 
0884   c1->Divide(4, 2);
0885   int idx = 1;
0886   TPad * p;
0887 
0888   for (int i = 0; i < beam_mom.size(); ++i)
0889     {
0890       p = (TPad *) c1->cd(idx++);
0891       c1->Update();
0892 
0893       const double momemtum = beam_mom[i];
0894       const TString histname = Form("hLineShape%.0fGeV", momemtum);
0895 
0896       TH1F * h = new TH1F(histname, histname + ";Observed energy (GeV)",
0897           (momemtum <= 6 ? 25 : 50), momemtum / 2, momemtum * 1.5);
0898       T->Draw(cluster_name + ".sum_E>>" + histname,
0899           Form("abs(abs(info.beam_mom)-%f)/%f<.06", momemtum, momemtum));
0900 
0901       h->Sumw2();
0902 
0903       TF1 * fgaus_g = new TF1("fgaus_LG_g", "gaus",
0904           h->GetMean() - 1 * h->GetRMS(), h->GetMean() + 4 * h->GetRMS());
0905       fgaus_g->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
0906           h->GetMean() + 2 * h->GetRMS());
0907       h->Fit(fgaus_g, "MR0N");
0908 
0909       TF1 * fgaus = new TF1("fgaus_LG", "gaus",
0910           fgaus_g->GetParameter(1) - 2 * fgaus_g->GetParameter(2),
0911           fgaus_g->GetParameter(1) + 2 * fgaus_g->GetParameter(2));
0912       fgaus->SetParameters(fgaus_g->GetParameter(0), fgaus_g->GetParameter(1),
0913           fgaus_g->GetParameter(2));
0914       h->Fit(fgaus, "MR");
0915 
0916       h->SetLineWidth(2);
0917       h->SetMarkerStyle(kFullCircle);
0918       fgaus->SetLineWidth(2);
0919 
0920       h->SetTitle(
0921           Form("%.0f GeV/c: #DeltaE/<E> = %.1f%%", momemtum,
0922               100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));
0923 
0924       mean.push_back(fgaus->GetParameter(1));
0925       mean_err.push_back(fgaus->GetParError(1));
0926       res.push_back(fgaus->GetParameter(2) / fgaus->GetParameter(1));
0927       res_err.push_back(fgaus->GetParError(2) / fgaus->GetParameter(1));
0928 
0929     }
0930 
0931   SaveCanvas(c1,
0932       TString(_file0->GetName()) + "_DrawPrototype2ShowerCalib_"
0933           + TString(c1->GetName()), kTRUE);
0934 
0935   TGraphErrors * ge_linear = new TGraphErrors(beam_mom.size(), &beam_mom[0],
0936       &mean[0], 0, &mean_err[0]);
0937 
0938   TGraphErrors * ge_res = new TGraphErrors(beam_mom.size(), &beam_mom[0],
0939       &res[0], 0, &res_err[0]);
0940   ge_res->GetHistogram()->SetStats(0);
0941   ge_res->Print();
0942 
0943   lin_res ret;
0944   ret.name = cluster_name;
0945   ret.linearity = ge_linear;
0946   ret.resolution = ge_res;
0947   ret.f_res = new TF1("f_calo_r", "sqrt([0]*[0]+[1]*[1]/x)/100", 0.5, 25);
0948   ge_res->Fit(ret.f_res, "RM0QN");
0949 
0950   ge_linear->SetLineColor(col);
0951   ge_linear->SetMarkerColor(col);
0952   ge_linear->SetLineWidth(2);
0953   ge_linear->SetMarkerStyle(kFullCircle);
0954   ge_linear->SetMarkerSize(1.5);
0955 
0956   ge_res->SetLineColor(col);
0957   ge_res->SetMarkerColor(col);
0958   ge_res->SetLineWidth(2);
0959   ge_res->SetMarkerStyle(kFullCircle);
0960   ge_res->SetMarkerSize(1.5);
0961 
0962   ge_res->GetHistogram()->SetStats(0);
0963 
0964   ret.f_res->SetLineColor(col);
0965   ret.f_res->SetLineWidth(3);
0966 
0967   return ret;
0968 }
0969 
0970 TGraphErrors *
0971 FitProfile(const TH2 * h2)
0972 {
0973 
0974   TProfile * p2 = h2->ProfileX();
0975 
0976   int n = 0;
0977   double x[1000];
0978   double ex[1000];
0979   double y[1000];
0980   double ey[1000];
0981 
0982   for (int i = 1; i <= h2->GetNbinsX(); i++)
0983     {
0984       TH1D * h1 = h2->ProjectionY(Form("htmp_%d", rand()), i, i);
0985 
0986       if (h1->GetSum() < 30)
0987         {
0988           cout << "FitProfile - ignore bin " << i << endl;
0989           continue;
0990         }
0991       else
0992         {
0993           cout << "FitProfile - fit bin " << i << endl;
0994         }
0995 
0996       TF1 fgaus("fgaus", "gaus", -p2->GetBinError(i) * 4,
0997           p2->GetBinError(i) * 4);
0998 
0999       TF1 f2(Form("dgaus"), "gaus + [3]*exp(-0.5*((x-[1])/[4])**2) + [5]",
1000           -p2->GetBinError(i) * 4, p2->GetBinError(i) * 4);
1001 
1002       fgaus.SetParameter(1, p2->GetBinContent(i));
1003       fgaus.SetParameter(2, p2->GetBinError(i));
1004 
1005       h1->Fit(&fgaus, "MQ");
1006 
1007       f2.SetParameters(fgaus.GetParameter(0) / 2, fgaus.GetParameter(1),
1008           fgaus.GetParameter(2), fgaus.GetParameter(0) / 2,
1009           fgaus.GetParameter(2) / 4, 0);
1010 
1011       h1->Fit(&f2, "MQ");
1012 
1013 //      new TCanvas;
1014 //      h1->Draw();
1015 //      fgaus.Draw("same");
1016 //      break;
1017 
1018       x[n] = p2->GetBinCenter(i);
1019       ex[n] = (p2->GetBinCenter(2) - p2->GetBinCenter(1)) / 2;
1020       y[n] = fgaus.GetParameter(1);
1021       ey[n] = fgaus.GetParError(1);
1022 
1023 //      p2->SetBinContent(i, fgaus.GetParameter(1));
1024 //      p2->SetBinError(i, fgaus.GetParameter(2));
1025 
1026       n++;
1027       delete h1;
1028     }
1029 
1030   return new TGraphErrors(n, x, y, ex, ey);
1031 }
1032
[ 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!