-
Notifications
You must be signed in to change notification settings - Fork 0
/
shapeToy.C
97 lines (83 loc) · 2.56 KB
/
shapeToy.C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
{
int nEvents=1000;
float lightYield=1000.;
float lightYieldResolution=0.1;
float tMean=50;
float tRes=4;
float triggerThreshold=40;
float decayTime=45;
//average current profile from PMT
TProfile full("full","full",1200,0,600);
TProfile fullTriggered("fullTriggered","fullTriggered",1200,0,600);
//scintillation decay function
TF1 *f=new TF1("f","TMath::Exp(-x/[0])",0,600);
f->SetParameter(0,decayTime);
TFile *fOut=new TFile("shapeToy.root","RECREATE");
fOut->mkdir("events");
fOut->cd("events");
for (int i=0;i<nEvents;++i)
{
int nph=gRandom->Gaus(lightYield,lightYieldResolution*lightYield);
//photo current for this event from the PMT (convoluting exp decay with gaussian PMT response)x
TH1F current("event","event",1200,0,600);
for (int ipho=0;ipho<nph;++ipho)
{
float t=f->GetRandom();
TF1 g("g","TMath::Gaus(x,[0],[1])",0,600);
g.SetParameter(0,tMean+t);
g.SetParameter(1,tRes);
for (int ib=1;ib<current.GetNbinsX();++ib)
{
float i=g.Eval(current.GetBinCenter(ib));
current.Fill(current.GetBinCenter(ib),i);
}
}
for (int ib=1;ib<current.GetNbinsX();++ib)
current.SetBinError(ib,0);
//find trigger time
TH1F* currentTriggered=(TH1F*)current.Clone("eventTriggered");
int triggerSample=-1;
for (int ib=1;ib<current.GetNbinsX();++ib)
if (current.GetBinContent(ib)>triggerThreshold)
{
triggerSample=ib-1;
break;
}
//skio events without a valid trigger
if (triggerSample==-1)
continue;
//optimised to have ~maximum of triggered function not too far off
triggerSample=triggerSample-60;
//shift current according to trigger time
for (int ib=1;ib<current.GetNbinsX();++ib)
{
if (ib+triggerSample<=current.GetNbinsX())
{
currentTriggered->SetBinContent(ib,currentTriggered->GetBinContent(ib+triggerSample));
currentTriggered->SetBinError(ib,0);
}
else
{
currentTriggered->SetBinContent(ib,0);
currentTriggered->SetBinError(ib,0);
}
}
//now fill the average function
for (int ib=1;ib<=current.GetNbinsX();++ib)
{
full.Fill(current.GetBinCenter(ib),current.GetBinContent(ib));
fullTriggered.Fill(currentTriggered->GetBinCenter(ib),currentTriggered->GetBinContent(ib));
}
//Save first 100 events
if (i<100)
{
current.Write(Form("event_%d",i));
currentTriggered->Write(Form("eventTriggered_%d",i));
}
delete currentTriggered;
}
fOut->cd("/");
full.Write();
fullTriggered.Write();
fOut->Close();
}