ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/HFmon2012/macros/simulateRaddam.C
Revision: 1.1
Committed: Mon Aug 13 14:57:43 2012 UTC (12 years, 8 months ago) by yilmaz
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Log Message:
update analysis

File Contents

# User Rev Content
1 yilmaz 1.1 #include <string>
2     #include <vector>
3     #include <iostream>
4     #include <map>
5     #include <utility>
6    
7     #include "HistogramManager.h"
8     #include "Map.h"
9    
10     #include "TNtuple.h"
11     #include "TRandom.h"
12    
13     #include "TTree.h"
14     #include "TFile.h"
15     #include "TH1I.h"
16     #include "TH1F.h"
17     #include "TH2F.h"
18     #include "TProfile.h"
19     #include "TString.h"
20     #include "TMath.h"
21     #include "TF1.h"
22     #include "TString.h"
23     #include "TCanvas.h"
24     #include "TCut.h"
25    
26     #include <fstream>
27     #include <vector>
28     #include <iomanip>
29     #include <iostream>
30    
31     #include "Laser.h"
32    
33     using namespace std;
34    
35    
36     TNtuple* simulate(double inputR = 0.5, double inputW = 0.1, double inputNoise = 10, const char* ntname = "nt50"){
37    
38     Laser *laz = new Laser(inputR,inputW, inputNoise);
39    
40    
41     TF1* f0 = new TF1("f0","pol0",-0.5,2.5);
42     f0->SetParameter(0,1);
43     // f0->SetParameter(1,4);
44     // f0->SetParameter(2,1);
45    
46     TNtuple* nt = new TNtuple(Form(ntname),"","s0:s1:s2:s3:sts0:sts1:sts2:sts3:sts4:sts5:sts6:sts7:sts8:sts9:ts:w:wts:r:rts:p");
47    
48     TH1D* hit = new TH1D("hit","",10,0,10);
49    
50     int Nev = 5000;
51    
52     for(int iev = 0; iev < Nev; iev++){
53    
54     hit->Reset();
55    
56     double p = f0->GetRandom();
57     laz->fillHit(hit,p);
58    
59     double s1 = 0., s2 = 0., s0 = 0, s3 = 0;
60     int its = 0;
61     findPeak(hit,s0,s1,s2,s3,its);
62    
63     double sts0 = hit->GetBinContent(1);
64     double sts1 = hit->GetBinContent(2);
65     double sts2 = hit->GetBinContent(3);
66     double sts3 = hit->GetBinContent(4);
67     double sts4 = hit->GetBinContent(5);
68     double sts5 = hit->GetBinContent(6);
69    
70     double sts6 = hit->GetBinContent(7);
71     double sts7 = hit->GetBinContent(8);
72     double sts8 = hit->GetBinContent(9);
73     double sts9 = hit->GetBinContent(10);
74    
75     double rts = sts4 / sts3;
76    
77     double wing = (s1+s2)/hit->Integral();
78     double wingts = (sts3+sts4)/hit->Integral();
79    
80     double r = s2/s1;
81    
82     float entry[100] = {s0,s1,s2,s3,sts0,sts1,sts2,sts3,sts4,sts5,sts6,sts7,sts8,sts9,its,wing,wingts,r,rts,p};
83     nt->Fill(entry);
84    
85    
86     }
87    
88    
89     // nt->Draw("r:p","r > -0.2 &&r < 1.2","");
90    
91     return nt;
92    
93     // hit->Draw();
94    
95     }
96    
97    
98     void simulateRaddam(){
99    
100     bool remake = 1;
101    
102     double r[100] = { 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.97, 0.98, 1. };
103     double w[100] = { 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5 };
104     double n[100] = { 0.5, 1, 2, 5, 10, 20, 100 };
105    
106     int Nr = 100;
107     int Nw = 10;
108     int Nn = 7;
109    
110     for(int i = 0; i < Nr; ++i){
111     r[i] = 1./Nr*(i+1);
112     }
113    
114    
115    
116     int Nplot = 3;
117    
118     TFile* outf, *inf;
119    
120     TNtuple* nt[20][20][20];
121     TH2D* h[20][20][20];
122     TH2D* hw[20][20][20];
123     TProfile* p[20][20][20];
124     TProfile* pw[20][20][20];
125    
126     TNtuple* ntx[100];
127     TH2D* hx[100];
128     TProfile* px[100];
129    
130     int id = 0;
131    
132     int color[100] = {1,2,3,4,5,6,7,8,9,10,11,12};
133    
134     TCut wing("w > 0.95");
135    
136     for(int i = 7; i < Nr; ++i){
137     for(int j = 0; j < Nw; ++j){
138     for(int k = 0; k < Nn; ++k){
139     if(remake) outf = new TFile(Form("simulations_%d_%d_%d.root",i,j,k),"recreate");
140    
141     if(remake){
142     nt[i][j][k] = simulate(r[i],w[j],n[k],Form("nt_%d_%d_%d",i,j,k));
143     }else{
144     nt[i][j][k] = (TNtuple*)inf->Get(Form("nt_%d_%d_%d",i,j,k));
145     }
146    
147     hw[i][j][k] = new TH2D(Form("hw_%d_%d_%d",i,j,k),"",200,-2,6,200,-0.2,1.2);
148     pw[i][j][k] = new TProfile(Form("pw_%d_%d_%d",i,j,k),"",200,-2,6);
149     nt[i][j][k]->Draw(Form("rts:p>>%s",hw[i][j][k]->GetName()),wing,"");
150     nt[i][j][k]->Draw(Form("rts:p>>%s",pw[i][j][k]->GetName()),wing&&"rts < 1.5 && rts > -0.5","");
151    
152     h[i][j][k] = new TH2D(Form("h_%d_%d_%d",i,j,k),"",200,-2,6,200,-0.2,1.2);
153     p[i][j][k] = new TProfile(Form("p_%d_%d_%d",i,j,k),"",200,-2,6);
154     nt[i][j][k]->Draw(Form("rts:p>>%s",h[i][j][k]->GetName()),"");
155     nt[i][j][k]->Draw(Form("rts:p>>%s",p[i][j][k]->GetName()),"rts < 1.5 && rts > -0.5","");
156    
157     // nt[i][j][k]->Delete();
158     outf->Write();
159     outf->Close();
160     id++;
161    
162     }
163     }
164     }
165    
166    
167     }
168    
169