| 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 |
|
