EMCal/EMCALmacros/ReadHis2.C

0001 #include <cmath>
0002 #include <TFile.h>
0003 #include <TString.h>
0004 #include <TLine.h>
0005 #include <TTree.h>
0006 #include <TLatex.h>
0007 #include <TGraphErrors.h>
0008 #include <cassert>
0009 #include "SaveCanvas.C"
0010 #include "SetOKStyle.C"
0011 #include <iostream>
0012 #include <fstream>
0013 #include "TROOT.h"
0014 #include "TH1.h"
0015 #include "TTree.h"
0016 using namespace std;
0017
0018
0019
0020 using std::cout;
0021 using std::endl;
0022 #endif
0023
0024
0025
0026
0027 void ReadHis2()
0028 {
0029
0030 int Method = 0;
0031
0032
0033   gSystem->Load("libg4eval.so");
0034   gSystem->Load("libqa_modules.so");
0035   gSystem->Load("libPrototype3.so");
0036  gStyle->SetOptFit(0);
0037    gStyle->SetOptStat(0);
0038
0039    double step = 1;
0040
0041 // Fitting: Method = 1; Mean: Method = 0//
0042
0043
0044    int Ini = 3543;
0045 int Final = 3579;
0046
0047
0048 double Xmin = 170;
0049 double Xmax = 310.0;
0050
0051 double Ymin = 60;
0052 double Ymax = 245.0;
0053
0054
0055 /*
0056 int Xmin = 170;
0057 int Xmax = 340.0;
0058
0059 int Ymin = 60;
0060 int Ymax = 245.0;
0061 */
0062
0063 double Energy;
0064 double Energy_err;
0065 double Width;
0066 double Width_err;
0067 double Center;
0068 double Center_err;
0069
0070 double x;
0071 double y;
0072 double xhigh;
0073 double yhigh;
0074 int binxhigh;
0075 int binyhigh;
0076
0077 int binx;
0078 int biny;
0079 int XBins = (Xmax - Xmin)/step;
0080 int YBins = (Ymax - Ymin)/step;
0081
0082 //int XBins = 200;
0083 //int YBins = 200;
0084 char Name[512];
0085
0086 int index = 0;
0087 double average;
0088 int indexmax = XBins * YBins;
0089 char Filename[512];
0090 char inputfile[512];
0091
0092 //sprintf(inputfile,"Hisfiles130/HisAll.root");
0093
0094 sprintf(inputfile,"His2.root");
0095 //sprintf(inputfile,"His3.root");
0096
0097 TFile *fin = new TFile(inputfile);
0098 TH2D *EnPo= new TH2D("EnPo","",XBins,Xmin,Xmax,YBins,Ymin,Ymax);
0099
0100 TF1 *f1 = new TF1("f1","gaus",6,11);
0101 cout << "XBins = " << XBins << endl;
0102 cout << "YBins = " << YBins << endl;
0103
0104 char hisname[512];
0105
0106 for(int i = 0; i < XBins; i++)
0107 {
0108
0109 x = Xmin + step*i;
0110
0111 xhigh = Xmin + step*(i+1);
0112
0113
0114 for(int j = 0; j < YBins; j++)
0115 {
0116
0117 TCanvas *c1 = new TCanvas("c1", "c1",0,0,800,600);
0118     y = Ymin + step*j;
0119 yhigh = Ymin + step*(j+1);
0120
0121 //etTitle("Total Energy vs Position by Fitting");
0122
0123
0124 binx = Energyhis->GetXaxis()->FindBin(x+0.0001);
0125 biny = Energyhis->GetYaxis()->FindBin(y+0.0001);
0126
0127 //binxhigh = binx +1;
0128 //binyhigh = biny +1;
0129
0130
0131
0132 binxhigh = Energyhis->GetXaxis()->FindBin(xhigh-0.0001);
0133 binyhigh = Energyhis->GetYaxis()->FindBin(yhigh-0.0001);
0134
0135
0136 TH1D *h5= new TH1D("h5","",25,4,10);
0137
0138
0139
0140
0141 h5 = Energyhis->ProjectionZ("h5",binx,binxhigh,biny,binyhigh);
0142
0143 //return;
0144
0145
0146
0147 //return;
0148 //
0149
0150 //Fitting Method //
0151 if(Method ==1){
0152
0153 //h5->Rebin(4);
0154
0155
0156 average = h5->GetMean();
0157
0158 if(average > 1)
0159 {
0160
0161  h5->Fit(f1,"R");
0162
0163 c1->Update();
0164
0165
0166
0167 Energy = f1->GetParameter(0);
0168 Energy_err = f1->GetParError(0);
0169 Center = f1->GetParameter(1);
0170 Center_err = f1->GetParError(1);
0171 Width = f1->GetParameter(2);
0172 Width_err = f1->GetParError(2);
0173
0174
0175 // Center = h5->GetMean();
0176
0177 // cout << "x = " <<  x << endl;
0178
0179 //cout << "y = " <<  y << endl;
0180
0181 if(Center > 11 || Center < 4){
0182
0183     cout << "Old Center = " << Center << endl;
0184
0185
0186 Center = h5->GetMean();
0187
0188 cout << "New Center = " << Center << endl;
0189 }
0190
0191
0192
0193 if(5.5 < Center < 6.5){
0194 h5->Draw();
0195
0196  h5->Fit(f1,"R");
0197
0198 sprintf(Name,"EMCAL -10 Degree Energy Distribution in x = %d mm and y = %d mm ",x,y)
0199 h5->GetXaxis()->SetTitle("Energy (GeV)");
0200 h5->GetYaxis()->SetTitle("Counts");
0201 h5->SetTitle(Name);
0202 c1->Update();
0203
0204
0205 sprintf(hisname,"pngfiles/Plot%d-%d.png",x,y);
0206
0207   c1->SaveAs(hisname);
0208 }
0209
0210 }
0211
0212
0213 if(average < 1) Center = 0;
0214
0215
0216 }
0217
0218
0219
0220 //Direct Mean Method//
0221
0222 if(Method == 0) Center = h5->GetMean();
0223
0224
0225 cout << "Mean =  "  << Center << endl;
0226
0227 EnPo->SetBinContent(binx,biny,Center);
0228
0229
0230 index = index + 1;
0231
0232 }
0233
0234 }
0235
0236 TCanvas *c22 = new TCanvas("c22", "c22",0,0,800,600);
0237
0238 c22->cd();
0239
0240 EnPo->Draw();
0241
0242     c22->Update();
0243
0244   c22->SaveAs("EnPo.png");
0245
0246
0247 //sprintf(Filename,"EnergyPosition.root");
0248
0249 sprintf(Filename,"EnergyPosition2.root");
0250 //sprintf(Filename,"EnergyPosition3.root");
0251
0252   TFile *fout = new TFile(Filename,"RECREATE");
0253 EnPo->Write();
0254
0255
0256 }