AlcapDAQ  1
EvdE.cpp
Go to the documentation of this file.
1 #include "EvdE.h"
2 #include "TAnalysedPulse.h"
3 #include "definitions.h"
4 
5 #include <iostream>
6 #include <cmath>
7 #include <string>
8 #include <vector>
9 #include <map>
10 
11 #include "TH1D.h"
12 #include "TH2D.h"
13 
14 extern StringAnalPulseMap gAnalysedPulseMap;
15 static int nSec; // Number of sections in the silicon thin detectors
16 static TH2D* hEvdE[2]; // Histograms (right, left)
17 static TH1D* hEvdE_log[2];
18 static TH1D* hEprotons[2];
19 static TH2D* hEvdEprotons[2];
20 
21 EvdE::EvdE(char *HistogramDirectoryName):FillHistBase(HistogramDirectoryName)
22 {;}
23 
24 EvdE::EvdE(char *HistogramDirectoryName, double t0, double t1) :
25  FillHistBase(HistogramDirectoryName)
26 {
27  // Thin silicon is in quadrants
28  nSec = 4;
29  // Time in ns
30  tCoincidence = 1000.;
31  tPP = 10000.;
32  tStart = t0;
33  tStop = t1;
34  // Energy cuts in ADC
35  adcCutThick[0] = 80.;
36  adcCutThick[1] = 80.;
37  adcCutThin[0] = 160.;
38  adcCutThin[1] = 160.;
39  eCutThin[0] = adcCutThin[0] * 2.6;
40  eCutThin[1] = adcCutThin[1] * 2.6;
41 
42  slopeThick[0] = 8.06;
43  offsetThick[0] = -192.70;
44  slopeThick[1] = 7.93;
45  offsetThick[1] = -116.33;
46 
47  slopeThin[0][0] = 2.57;
48  offsetThin[0][0] = -87.96;
49  slopeThin[0][1] = 2.72;
50  offsetThin[0][1] = -292.15;
51  slopeThin[0][2] = 2.52;
52  offsetThin[0][2] = -30.19;
53  slopeThin[0][3] = 2.57;
54  offsetThin[0][3] = -103.67;
55 
56  slopeThin[1][0] = 2.66;
57  offsetThin[1][0] = -117.34;
58  slopeThin[1][1] = 2.56;
59  offsetThin[1][1] = -44.45;
60  slopeThin[1][2] = 2.68;
61  offsetThin[1][2] = -129.73;
62  slopeThin[1][3] = 2.58;
63  offsetThin[1][3] = -85.48;
64 
65  // Prepare histograms
66  hEvdE[0] = new TH2D("hEvdE_Right", "dEdx vs E (Right);E + dE (keV);dE (keV)",
67  1000, 0., 14000.,
68  500, 0., 7000.);
69  hEvdE[1] = new TH2D("hEvdE_Left", "dEdx vs E (Left);E + dE (keV);dE (keV)",
70  1000, 0., 14000.,
71  500, 0., 7000.);
72 
73  hEvdE_log[0] = new TH1D("hEvdE_Right_log", "logdE +logE (Right)", 100, 10., 30.);
74  hEvdE_log[1] = new TH1D("hEvdE_Left_log", "logdE +logE (Left)", 100, 10., 30.);
75 
76  hEprotons[1] = new TH1D("hEp_Left",
77  "Energy of protons (Left)", 1000, 0., 14000.);
78  hEprotons[0] = new TH1D("hEp_Right",
79  "Energy of protons (Right)", 1000, 0., 14000.);
80 
81  hEvdEprotons[0] = new TH2D("hEvdE_Proton_Right", "dEdx vs E (Right);E + dE (keV);dE (keV)",
82  1000, 0., 14000.,
83  500, 0., 7000.);
84  hEvdEprotons[1] = new TH2D("hEvdE_Proton_Left", "dEdx vs E (Left);E + dE (keV);dE (keV)",
85  1000, 0., 14000.,
86  500, 0., 7000.);
87 
88  gDirectory->cd("/");
89 }
90 
91 EvdE::~EvdE(){
92 }
93 
94 int EvdE::ProcessEntry(TGlobalData *gData, TSetupData *gSetup) {
95 
96  // Look for muSc and at least one arm
97  static AnalysedPulseList *musc;
98  static AnalysedPulseList *sirThin[4], *silThin[4];
99  static AnalysedPulseList *sirThick, *silThick;
100  if (gAnalysedPulseMap.count("muSc")) musc = &gAnalysedPulseMap.at("muSc");
101  else return 0;
102 
103  if (gAnalysedPulseMap.count("SiR1-1-S")) sirThin[0] = &gAnalysedPulseMap.at("SiR1-1-S");
104  else sirThin[0] = NULL;
105  if (gAnalysedPulseMap.count("SiR1-2-S")) sirThin[1] = &gAnalysedPulseMap.at("SiR1-2-S");
106  else sirThin[1] = NULL;
107  if (gAnalysedPulseMap.count("SiR1-3-S")) sirThin[2] = &gAnalysedPulseMap.at("SiR1-3-S");
108  else sirThin[2] = NULL;
109  if (gAnalysedPulseMap.count("SiR1-4-S")) sirThin[3] = &gAnalysedPulseMap.at("SiR1-4-S");
110  else sirThin[3] = NULL;
111 
112  if (gAnalysedPulseMap.count("SiL1-1-S")) silThin[0] = &gAnalysedPulseMap.at("SiL1-1-S");
113  else silThin[0] = NULL;
114  if (gAnalysedPulseMap.count("SiL1-2-S")) silThin[1] = &gAnalysedPulseMap.at("SiL1-2-S");
115  else silThin[1] = NULL;
116  if (gAnalysedPulseMap.count("SiL1-3-S")) silThin[2] = &gAnalysedPulseMap.at("SiL1-3-S");
117  else silThin[2] = NULL;
118  if (gAnalysedPulseMap.count("SiL1-4-S")) silThin[3] = &gAnalysedPulseMap.at("SiL1-4-S");
119  else silThin[3] = NULL;
120 
121  if (gAnalysedPulseMap.count("SiR2-S")) sirThick = &gAnalysedPulseMap.at("SiR2-S");
122  else sirThick = NULL;
123 
124  if (gAnalysedPulseMap.count("SiL2-S")) silThick = &gAnalysedPulseMap.at("SiL2-S");
125  else silThick = NULL;
126 
127  // If there's nothing logical to get data from, don't bother analysing
128  if (!(sirThick && (sirThin[0] || sirThin[1] || sirThin[2] || sirThin[3])) &&
129  !(silThick && (silThin[0] || silThin[1] || silThin[2] || silThin[3])))
130  return 0;
131 
132 
133  // Iterators through aforementioned vectors
134  static int iSec;
135  static AnalysedPulseList::iterator cMusc;
136  static AnalysedPulseList::iterator cSiRThick, cSiLThick;
137  static AnalysedPulseList::iterator cSiRThin[4], cSiLThin[4];
138  cMusc = musc->begin();
139  if (sirThick)
140  cSiRThick = sirThick->begin();
141  if (silThick)
142  cSiLThick = silThick->begin();
143  for (iSec = 0; iSec < nSec; ++iSec)
144  {
145  if (sirThin[iSec])
146  cSiRThin[iSec] = (sirThin[iSec])->begin();
147 
148  if (silThin[iSec])
149  cSiLThin[iSec] = (silThin[iSec])->begin();
150  }
151 
152  static double tMusc;
153  static std::vector<double> tSiRThick, tSiLThick;
154  static std::vector<double> eSiRThick, eSiLThick;
155  static unsigned int iHit, nThickHits;
156 
157  for (; cMusc != musc->end(); ++cMusc)
158  {
159  tMusc = (*cMusc)->GetTime();
160  tSiRThick.clear(); tSiLThick.clear();
161  eSiRThick.clear(); eSiLThick.clear();
162  // If there's pile-up, advance beyond it
163  if ((cMusc+1) != musc->end() && ((*(cMusc + 1))->GetTime() - tMusc) < tPP)
164  {
165  ++cMusc;
166  continue;
167  }
168 
169  // Advance silicon pulses beyond time of muSc
170  if (sirThick)
171  while (cSiRThick != sirThick->end()
172  && ((*cSiRThick)->GetTime() < (tMusc + tStart)))
173  ++cSiRThick;
174 
175  if (silThick)
176  while ((cSiLThick != silThick->end())
177  && ((*cSiLThick)->GetTime() < tMusc + tStart))
178  ++cSiLThick;
179 
180  for (iSec = 0; iSec < nSec; ++iSec)
181  {
182  if (sirThin[iSec])
183  while (cSiRThin[iSec] != sirThin[iSec]->end()
184  && ((*(cSiRThin[iSec]))->GetTime() < tMusc + tStart))
185  ++(cSiRThin[iSec]);
186 
187  if (silThin[iSec])
188  while (cSiLThin[iSec] != silThin[iSec]->end()
189  && ((*(cSiLThin[iSec]))->GetTime() < tMusc + tStart))
190  ++(cSiLThin[iSec]);
191  }
192 
193  // Check for hits
194  // Look for hits in thick detectors first
195  if (sirThick) // If there was at least one hit in the thick detector
196  {
197  while (cSiRThick != sirThick->end()
198  && (((*cSiRThick)->GetTime() - tMusc) < tStop)) // Look for hits in the decay window after muSc hit
199  {
200  if ((*cSiRThick)->GetEnergy() > adcCutThick[0])
201  {// If a hit passes the amplitude cut
202  tSiRThick.push_back((*cSiRThick)->GetTime()); // Save the time
203  eSiRThick.push_back((*cSiRThick)->GetEnergy()); // And the amplitude
204  }
205  ++cSiRThick;
206  }
207  }
208 
209  if (silThick)
210  {
211  while (cSiLThick != silThick->end()
212  && (((*cSiLThick)->GetTime() - tMusc < tStop)))
213  {
214  if ((*cSiLThick)->GetEnergy() > adcCutThick[1])
215  {
216  tSiLThick.push_back((*cSiLThick)->GetTime());
217  eSiLThick.push_back((*cSiLThick)->GetEnergy());
218  }
219  ++cSiLThick;
220  }
221  }
222 
223  if (tSiRThick.size() > 0) // If there was at least one hit in the thick window
224  {
225  nThickHits = tSiRThick.size();
226  for (iHit = 0; iHit < nThickHits; ++iHit)
227  {// Loop through the times of all of the hits
228  double dE = 0.; // Keep track of energy
229  for (iSec = 0; iSec < nSec; ++iSec)// Loop through all the sections of the thin silicon
230  {
231  if (sirThin[iSec])
232  {// And, if there are hits in this thin detector
233  while ((cSiRThin[iSec] != sirThin[iSec]->end()) &&
234  (*(cSiRThin[iSec]))->GetTime() < (tSiRThick[iHit] - tCoincidence)) // Catch the iterator up to the beginning of the time window
235  ++(cSiRThin[iSec]);
236 
237  while (cSiRThin[iSec] != sirThin[iSec]->end() &&
238  (*(cSiRThin[iSec]))->GetTime() < tSiRThick[iHit] + tCoincidence)
239  { // And look for hits in this window
240  dE += offsetThin[0][iSec] + // Increment the energy
241  ((*(cSiRThin[iSec]))->GetEnergy())*slopeThin[0][iSec];
242  ++(cSiRThin[iSec]);
243  }
244  }
245  }
246  if (dE> eCutThin[0]) // If the energy is greater than our cut
247  {
248  double E = eSiRThick[iHit]*slopeThick[0] + offsetThick[0] + dE;
249  hEvdE[0]->Fill(E, dE);// Record total energy deposited in thin detector
250  hEvdE_log[0]->Fill(log(E) + log(dE));
251  if ((log(E) + log(dE) >14.5) && (log(E) + log(dE)<15.6))
252  {
253  hEprotons[0]->Fill(E);
254  hEvdEprotons[0]->Fill(E, dE);
255  }
256  }
257  }
258  }
259 
260  if (tSiLThick.size() > 0)
261  {
262  nThickHits = tSiLThick.size();
263  for (iHit = 0; iHit < nThickHits; ++iHit)
264  {
265  double dE = 0.;
266  for (iSec = 0; iSec < nSec; ++iSec)
267  {
268  if (silThin[iSec]) {
269  while (cSiLThin[iSec] != silThin[iSec]->end() &&
270  (*(cSiLThin[iSec]))->GetTime() < (tSiLThick[iHit]-tCoincidence))
271  ++(cSiLThin[iSec]);
272 
273  while (cSiLThin[iSec] != silThin[iSec]->end() &&
274  (*(cSiLThin[iSec]))->GetTime() < tSiLThick[iHit] + tCoincidence)
275  {
276  dE += offsetThin[1][iSec] + // Increment the energy
277  ((*(cSiLThin[iSec]))->GetEnergy())*slopeThin[1][iSec];
278  ++(cSiLThin[iSec]);
279  }
280  }
281  }
282  if (dE > eCutThin[1])
283  {
284  double E = eSiLThick[iHit]*slopeThick[1] + offsetThick[1] + dE;
285  hEvdE[1]->Fill(E, dE);
286  hEvdE_log[1]->Fill(log(E) + log(dE));
287  if ((log(E) + log(dE) >14.5) && (log(E) + log(dE)<15.6))
288  {
289  hEprotons[1]->Fill(E);
290  hEvdEprotons[1]->Fill(E, dE);
291  }
292  }
293  }
294  }
295  }
296  return 0;
297 }
Generated by   doxygen 1.8.4