HFmon2012/macros/getR.C



void concentrate(TH2D*& hh){

  TH2D* hm = (TH2D*)hh->Clone(Form("%s_merged",hh->GetName()));
  hm->Reset();

  int Nx = hh->GetNbinsX();
  int Ny = hh->GetNbinsY();

  const int Nmerge = 10;

  TH2D* hp[Nmerge];

  for(int i = 0; i < Nmerge; ++i){
    hp[i] = (TH2D*)hm->Clone(Form("hb%d",i));
    
    for(int ix = 0; ix < Nx+1; ++ix){
      for(int iy = 0; iy < Ny+1; ++iy){
        
        int ibin = hp[i]->GetBin(ix,iy);
        int ibinshift = hp[i]->GetBin(ix+(i*25),iy);

        double y = hh->GetBinContent(ibinshift);
        double e = hh->GetBinError(ibinshift);

        hp[i]->SetBinContent(ibin,y);
        hp[i]->SetBinError(ibin,e);
        
      }
    }

  }


  for(int i = 0; i < Nmerge; ++i){
    hm->Add(hp[i]);

  }

  hh = hm;

}


double getRwithTime(int run = 178579, int ieta = 30, int iphi = 57, int idepth = 1){

  double rangeMin = 1052;

  if(0){
    run = 187937;
    ieta = 30;
    iphi = 21;
    idepth = 1;

  };


  double rangeMax = rangeMin + 24;

  TH1::SetDefaultSumw2();
  TH2::SetDefaultSumw2();

  double r = 0;

  TFile* inf = new TFile(Form("../test/res%d.root",run));
  TNtuple * nt = (TNtuple*)inf->Get("nt");

  TCut channel(Form("ieta==%d&&iphi==%d&&idepth==%d",ieta,iphi,idepth));
  TCut reasonable("ry > -0.2 && ry < 2");
  TCut wing("(s3+s0)<0.10*(s1+s2)");

  TCut cut = channel && wing && reasonable;

  TH2D* h2 = new TH2D(Form("h2_%d_%d_%d_%d",run,ieta,iphi,idepth),";phase;R",250,1050,1300,100,0.05,2.05);
  TProfile* h = 0; //new TProfile("h",";phase;R",250,1050,1300);
  
  nt->SetAlias("ry","s1/s2");
  nt->SetAlias("phase","p-5");

  nt->Draw(Form("ry:phase>>%s",h2->GetName()),cut,"colz");
  //  nt->Draw("ry:p>>h",cut,"");

  concentrate(h2);

  h = h2->ProfileX();

  int iMin = 0;

  h->SetAxisRange(rangeMin,rangeMax);
  iMin = h->GetMaximumBin();
  rangeMin = h->GetBinLowEdge(iMin);

  if(rangeMin > 1070) rangeMin -= 25;
  rangeMax = rangeMin + 24;

  TF1* f = new TF1("f","[0]-[1]*TMath::Erf((x-[2])/[3])",rangeMin,rangeMax);
  f->SetParLimits(2,rangeMin,rangeMax);
  f->SetParLimits(3,5,80);
  f->SetParLimits(0,1,3);


  h->Fit(f,"R");

  h2->SetAxisRange(rangeMin,rangeMax);
  h2->Draw("colz");
  h->Draw("same");

  r = f->GetParameter(0) - f->GetParameter(1);

  cout<<"Result : "<<r<<endl;


  return r;


}


double getRwithoutTime(int run = 178579, int ieta = 30, int iphi = 57, int idepth = 1){

  TH1::SetDefaultSumw2();
  TH2::SetDefaultSumw2();

  double r = 0;

  TFile* inf = new TFile(Form("../test/res%d.root",run));
  TNtuple * nt = (TNtuple*)inf->Get("nt");

  TCut channel(Form("ieta==%d&&iphi==%d&&idepth==%d",ieta,iphi,idepth));
  TCut reasonable("ry > -0.2 && ry < 2");

  TCut wing3("(s3+s0)<0.03*(s1+s2)");
  TCut wing5("(s3+s0)<0.05*(s1+s2)");
  TCut wing10("(s3+s0)<0.10*(s1+s2)");

  TCut wing = wing3;

  TCut cut = channel && wing;

  TF1* f = new TF1(Form("fg_%d_%d_%d_%d",run,ieta,iphi,idepth),"gaus(0)+[3]",0,2);

  int Nbin = 50;
  TH1D* h = new TH1D(Form("hg_%d_%d_%d_%d",run,ieta,iphi,idepth),"",Nbin,0,2);
  TH1D* h3 = new TH1D(Form("hg3_%d_%d_%d_%d",run,ieta,iphi,idepth),"",Nbin,0,2);
  TH1D* h5 = new TH1D(Form("hg5_%d_%d_%d_%d",run,ieta,iphi,idepth),"",Nbin,0,2);
  TH1D* h10 = new TH1D(Form("hg10_%d_%d_%d_%d",run,ieta,iphi,idepth),"",Nbin,0,2);

  nt->SetAlias("ry","s1/s2");
  if(0){
    nt->Draw(Form("ry>>%s",h3->GetName()),cut&&wing3,"");
    nt->Draw(Form("ry>>%s",h5->GetName()),cut&&wing5,"");
    nt->Draw(Form("ry>>%s",h10->GetName()),cut&&wing10,"");
  }
  nt->Draw(Form("ry>>%s",h->GetName()),cut,"");

  f->SetParameter(0,h->GetMaximum());
  f->SetParameter(1,h->GetMean());
  f->SetParameter(2,h->GetRMS());


  h->Fit(f,"R");

  r = f->GetParameter(1);

  new TCanvas();  
  return r;

}


double getR(int run = 130158, int ieta = 30, int iphi = 57, int idepth = 1, bool useTime = 1){


  double r2 = getRwithoutTime(run,ieta,iphi,idepth);
  double r1 = getRwithTime(run,ieta,iphi,idepth);

  TLine* line = new TLine(1000,r2,1175,r2);  
  line->SetLineWidth(5);
  line->SetLineColor(1);
  line->SetLineStyle(2);

  line->Draw("same");
  if(useTime) return r1;
  else return r2;

}


Revision:	1.7
Committed:	Wed Aug 15 17:44:46 2012 UTC (12 years, 8 months ago) by makbiyik
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.6:	+33 -6 lines
Log Message:	super update
#	User	Rev	Content
1	yilmaz	1.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	makbiyik	1.3	double getRwithTime(int run = 178579, int ieta = 30, int iphi = 57, int idepth = 1){
49	yilmaz	1.1
50	makbiyik	1.3	double rangeMin = 1052;
51	yilmaz	1.1
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	makbiyik	1.3	TFile* inf = new TFile(Form("../test/res%d.root",run));
69	makbiyik	1.5	TNtuple * nt = (TNtuple*)inf->Get("nt");
70	yilmaz	1.1
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	makbiyik	1.7
96			if(rangeMin > 1070) rangeMin -= 25;
97	yilmaz	1.1	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	makbiyik	1.3
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	makbiyik	1.6	TFile* inf = new TFile(Form("../test/res%d.root",run));
130	makbiyik	1.3	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	makbiyik	1.7
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	makbiyik	1.3
141	makbiyik	1.6	TCut cut = channel && wing;
142
143	makbiyik	1.7	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	makbiyik	1.3
151			nt->SetAlias("ry","s1/s2");
152	makbiyik	1.7	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	makbiyik	1.6	nt->Draw(Form("ry>>%s",h->GetName()),cut,"");
158	makbiyik	1.3
159	makbiyik	1.7	f->SetParameter(0,h->GetMaximum());
160			f->SetParameter(1,h->GetMean());
161			f->SetParameter(2,h->GetRMS());
162
163
164	makbiyik	1.6	h->Fit(f,"R");
165	makbiyik	1.3
166	makbiyik	1.6	r = f->GetParameter(1);
167	makbiyik	1.7
168			new TCanvas();
169	makbiyik	1.6	return r;
170	makbiyik	1.3
171			}
172
173
174
175
176
177
178	makbiyik	1.7	double getR(int run = 130158, int ieta = 30, int iphi = 57, int idepth = 1, bool useTime = 1){
179	makbiyik	1.6
180	makbiyik	1.3
181	makbiyik	1.6	double r2 = getRwithoutTime(run,ieta,iphi,idepth);
182			double r1 = getRwithTime(run,ieta,iphi,idepth);
183	makbiyik	1.3
184	makbiyik	1.7	TLine* line = new TLine(1000,r2,1175,r2);
185			line->SetLineWidth(5);
186			line->SetLineColor(1);
187			line->SetLineStyle(2);
188
189	makbiyik	1.6	line->Draw("same");
190			if(useTime) return r1;
191			else return r2;
192	makbiyik	1.3
193			}
194
195