Back to home page

sPhenix code displayed by LXR
 
 

    

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

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