Back to home page

sPhenix code displayed by LXR
 
 

    

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

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