| 1 |
|
| 2 |
|
| 3 |
|
| 4 |
|
| 5 |
void concentrate(TH2D*& hh){
|
| 6 |
|
| 7 |
TH2D* hm = (TH2D*)hh->Clone(Form("%s_merged",hh->GetName()));
|
| 8 |
hm->Reset();
|
| 9 |
|
| 10 |
int Nx = hh->GetNbinsX();
|
| 11 |
int Ny = hh->GetNbinsY();
|
| 12 |
|
| 13 |
const int Nmerge = 10;
|
| 14 |
|
| 15 |
TH2D* hp[Nmerge];
|
| 16 |
|
| 17 |
for(int i = 0; i < Nmerge; ++i){
|
| 18 |
hp[i] = (TH2D*)hm->Clone(Form("hb%d",i));
|
| 19 |
|
| 20 |
for(int ix = 0; ix < Nx+1; ++ix){
|
| 21 |
for(int iy = 0; iy < Ny+1; ++iy){
|
| 22 |
|
| 23 |
int ibin = hp[i]->GetBin(ix,iy);
|
| 24 |
int ibinshift = hp[i]->GetBin(ix+(i*25),iy);
|
| 25 |
|
| 26 |
double y = hh->GetBinContent(ibinshift);
|
| 27 |
double e = hh->GetBinError(ibinshift);
|
| 28 |
|
| 29 |
hp[i]->SetBinContent(ibin,y);
|
| 30 |
hp[i]->SetBinError(ibin,e);
|
| 31 |
|
| 32 |
}
|
| 33 |
}
|
| 34 |
|
| 35 |
}
|
| 36 |
|
| 37 |
|
| 38 |
for(int i = 0; i < Nmerge; ++i){
|
| 39 |
hm->Add(hp[i]);
|
| 40 |
|
| 41 |
}
|
| 42 |
|
| 43 |
hh = hm;
|
| 44 |
|
| 45 |
}
|
| 46 |
|
| 47 |
|
| 48 |
double getRwithTime(int run = 178579, int ieta = 30, int iphi = 57, int idepth = 1){
|
| 49 |
|
| 50 |
double rangeMin = 1052;
|
| 51 |
|
| 52 |
if(0){
|
| 53 |
run = 187937;
|
| 54 |
ieta = 30;
|
| 55 |
iphi = 21;
|
| 56 |
idepth = 1;
|
| 57 |
|
| 58 |
};
|
| 59 |
|
| 60 |
|
| 61 |
double rangeMax = rangeMin + 24;
|
| 62 |
|
| 63 |
TH1::SetDefaultSumw2();
|
| 64 |
TH2::SetDefaultSumw2();
|
| 65 |
|
| 66 |
double r = 0;
|
| 67 |
|
| 68 |
TFile* inf = new TFile(Form("../test/res%d.root",run));
|
| 69 |
TNtuple * nt = (TNtuple*)inf->Get("nt");
|
| 70 |
|
| 71 |
TCut channel(Form("ieta==%d&&iphi==%d&&idepth==%d",ieta,iphi,idepth));
|
| 72 |
TCut reasonable("ry > -0.2 && ry < 2");
|
| 73 |
TCut wing("(s3+s0)<0.10*(s1+s2)");
|
| 74 |
|
| 75 |
TCut cut = channel && wing && reasonable;
|
| 76 |
|
| 77 |
TH2D* h2 = new TH2D(Form("h2_%d_%d_%d_%d",run,ieta,iphi,idepth),";phase;R",250,1050,1300,100,0.05,2.05);
|
| 78 |
TProfile* h = 0; //new TProfile("h",";phase;R",250,1050,1300);
|
| 79 |
|
| 80 |
nt->SetAlias("ry","s1/s2");
|
| 81 |
nt->SetAlias("phase","p-5");
|
| 82 |
|
| 83 |
nt->Draw(Form("ry:phase>>%s",h2->GetName()),cut,"colz");
|
| 84 |
// nt->Draw("ry:p>>h",cut,"");
|
| 85 |
|
| 86 |
concentrate(h2);
|
| 87 |
|
| 88 |
h = h2->ProfileX();
|
| 89 |
|
| 90 |
int iMin = 0;
|
| 91 |
|
| 92 |
h->SetAxisRange(rangeMin,rangeMax);
|
| 93 |
iMin = h->GetMaximumBin();
|
| 94 |
rangeMin = h->GetBinLowEdge(iMin);
|
| 95 |
|
| 96 |
if(rangeMin > 1070) rangeMin -= 25;
|
| 97 |
rangeMax = rangeMin + 24;
|
| 98 |
|
| 99 |
TF1* f = new TF1("f","[0]-[1]*TMath::Erf((x-[2])/[3])",rangeMin,rangeMax);
|
| 100 |
f->SetParLimits(2,rangeMin,rangeMax);
|
| 101 |
f->SetParLimits(3,5,80);
|
| 102 |
f->SetParLimits(0,1,3);
|
| 103 |
|
| 104 |
|
| 105 |
h->Fit(f,"R");
|
| 106 |
|
| 107 |
h2->SetAxisRange(rangeMin,rangeMax);
|
| 108 |
h2->Draw("colz");
|
| 109 |
h->Draw("same");
|
| 110 |
|
| 111 |
r = f->GetParameter(0) - f->GetParameter(1);
|
| 112 |
|
| 113 |
cout<<"Result : "<<r<<endl;
|
| 114 |
|
| 115 |
|
| 116 |
return r;
|
| 117 |
|
| 118 |
|
| 119 |
}
|
| 120 |
|
| 121 |
|
| 122 |
double getRwithoutTime(int run = 178579, int ieta = 30, int iphi = 57, int idepth = 1){
|
| 123 |
|
| 124 |
TH1::SetDefaultSumw2();
|
| 125 |
TH2::SetDefaultSumw2();
|
| 126 |
|
| 127 |
double r = 0;
|
| 128 |
|
| 129 |
TFile* inf = new TFile(Form("../test/res%d.root",run));
|
| 130 |
TNtuple * nt = (TNtuple*)inf->Get("nt");
|
| 131 |
|
| 132 |
TCut channel(Form("ieta==%d&&iphi==%d&&idepth==%d",ieta,iphi,idepth));
|
| 133 |
TCut reasonable("ry > -0.2 && ry < 2");
|
| 134 |
|
| 135 |
TCut wing3("(s3+s0)<0.03*(s1+s2)");
|
| 136 |
TCut wing5("(s3+s0)<0.05*(s1+s2)");
|
| 137 |
TCut wing10("(s3+s0)<0.10*(s1+s2)");
|
| 138 |
|
| 139 |
TCut wing = wing3;
|
| 140 |
|
| 141 |
TCut cut = channel && wing;
|
| 142 |
|
| 143 |
TF1* f = new TF1(Form("fg_%d_%d_%d_%d",run,ieta,iphi,idepth),"gaus(0)+[3]",0,2);
|
| 144 |
|
| 145 |
int Nbin = 50;
|
| 146 |
TH1D* h = new TH1D(Form("hg_%d_%d_%d_%d",run,ieta,iphi,idepth),"",Nbin,0,2);
|
| 147 |
TH1D* h3 = new TH1D(Form("hg3_%d_%d_%d_%d",run,ieta,iphi,idepth),"",Nbin,0,2);
|
| 148 |
TH1D* h5 = new TH1D(Form("hg5_%d_%d_%d_%d",run,ieta,iphi,idepth),"",Nbin,0,2);
|
| 149 |
TH1D* h10 = new TH1D(Form("hg10_%d_%d_%d_%d",run,ieta,iphi,idepth),"",Nbin,0,2);
|
| 150 |
|
| 151 |
nt->SetAlias("ry","s1/s2");
|
| 152 |
if(0){
|
| 153 |
nt->Draw(Form("ry>>%s",h3->GetName()),cut&&wing3,"");
|
| 154 |
nt->Draw(Form("ry>>%s",h5->GetName()),cut&&wing5,"");
|
| 155 |
nt->Draw(Form("ry>>%s",h10->GetName()),cut&&wing10,"");
|
| 156 |
}
|
| 157 |
nt->Draw(Form("ry>>%s",h->GetName()),cut,"");
|
| 158 |
|
| 159 |
f->SetParameter(0,h->GetMaximum());
|
| 160 |
f->SetParameter(1,h->GetMean());
|
| 161 |
f->SetParameter(2,h->GetRMS());
|
| 162 |
|
| 163 |
|
| 164 |
h->Fit(f,"R");
|
| 165 |
|
| 166 |
r = f->GetParameter(1);
|
| 167 |
|
| 168 |
new TCanvas();
|
| 169 |
return r;
|
| 170 |
|
| 171 |
}
|
| 172 |
|
| 173 |
|
| 174 |
|
| 175 |
|
| 176 |
|
| 177 |
|
| 178 |
double getR(int run = 130158, int ieta = 30, int iphi = 57, int idepth = 1, bool useTime = 1){
|
| 179 |
|
| 180 |
|
| 181 |
double r2 = getRwithoutTime(run,ieta,iphi,idepth);
|
| 182 |
double r1 = getRwithTime(run,ieta,iphi,idepth);
|
| 183 |
|
| 184 |
TLine* line = new TLine(1000,r2,1175,r2);
|
| 185 |
line->SetLineWidth(5);
|
| 186 |
line->SetLineColor(1);
|
| 187 |
line->SetLineStyle(2);
|
| 188 |
|
| 189 |
line->Draw("same");
|
| 190 |
if(useTime) return r1;
|
| 191 |
else return r2;
|
| 192 |
|
| 193 |
}
|
| 194 |
|
| 195 |
|
