1 | #include <TStyle.h>
|
---|
2 | #include <TCanvas.h>
|
---|
3 | #include <TSystem.h>
|
---|
4 | #include <TF1.h>
|
---|
5 | #include <TProfile.h>
|
---|
6 | #include <TProfile2D.h>
|
---|
7 | #include <TMath.h>
|
---|
8 | #include <TFitResultPtr.h>
|
---|
9 | #include <TFitResult.h>
|
---|
10 |
|
---|
11 | #include "MH.h"
|
---|
12 | #include "MLog.h"
|
---|
13 | #include "MLogManip.h"
|
---|
14 | #include "MDirIter.h"
|
---|
15 | #include "MFillH.h"
|
---|
16 | #include "MEvtLoop.h"
|
---|
17 | #include "MCamEvent.h"
|
---|
18 | #include "MGeomApply.h"
|
---|
19 | #include "MTaskList.h"
|
---|
20 | #include "MParList.h"
|
---|
21 | #include "MContinue.h"
|
---|
22 | #include "MBinning.h"
|
---|
23 | #include "MHDrsCalib.h"
|
---|
24 | #include "MDrsCalibApply.h"
|
---|
25 | #include "MRawFitsRead.h"
|
---|
26 | #include "MBadPixelsCam.h"
|
---|
27 | #include "MStatusDisplay.h"
|
---|
28 | #include "MTaskInteractive.h"
|
---|
29 | #include "MPedestalSubtractedEvt.h"
|
---|
30 | #include "MHCamera.h"
|
---|
31 | #include "MGeomCamFACT.h"
|
---|
32 | #include "MRawRunHeader.h"
|
---|
33 |
|
---|
34 | using namespace std;
|
---|
35 |
|
---|
36 | MPedestalSubtractedEvt *fEvent = 0;
|
---|
37 | MLog *fLog = &gLog;
|
---|
38 |
|
---|
39 | struct Single
|
---|
40 | {
|
---|
41 | float fSignal;
|
---|
42 | float fTime;
|
---|
43 | };
|
---|
44 |
|
---|
45 | class MSingles : public MParContainer, public MCamEvent
|
---|
46 | {
|
---|
47 | Int_t fIntegrationWindow;
|
---|
48 |
|
---|
49 | vector<vector<Single>> fData;
|
---|
50 |
|
---|
51 | public:
|
---|
52 | MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(0)
|
---|
53 | {
|
---|
54 | fName = name ? name : "MSingles";
|
---|
55 | }
|
---|
56 |
|
---|
57 | void InitSize(const UInt_t i)
|
---|
58 | {
|
---|
59 | fData.resize(i);
|
---|
60 | }
|
---|
61 |
|
---|
62 | vector<Single> &operator[](UInt_t i) { return fData[i]; }
|
---|
63 | vector<Single> &GetVector(UInt_t i) { return fData[i]; }
|
---|
64 |
|
---|
65 | UInt_t GetNumPixels() const { return fData.size(); }
|
---|
66 |
|
---|
67 | void SetIntegrationWindow(Int_t w) { fIntegrationWindow = w; }
|
---|
68 | Int_t GetIntegrationWindow() const { return fIntegrationWindow; }
|
---|
69 |
|
---|
70 | Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const
|
---|
71 | {
|
---|
72 | return kTRUE;
|
---|
73 | }
|
---|
74 | void DrawPixelContent(Int_t num) const { }
|
---|
75 |
|
---|
76 | ClassDef(MSingles, 1)
|
---|
77 | };
|
---|
78 |
|
---|
79 | class MH2F : public TH2F
|
---|
80 | {
|
---|
81 | public:
|
---|
82 | Int_t Fill(Double_t x,Double_t y)
|
---|
83 | {
|
---|
84 | Int_t binx, biny, bin;
|
---|
85 | fEntries++;
|
---|
86 | binx = TMath::Nint(x+1);
|
---|
87 | biny = fYaxis.FindBin(y);
|
---|
88 | bin = biny*(fXaxis.GetNbins()+2) + binx;
|
---|
89 | AddBinContent(bin);
|
---|
90 | if (fSumw2.fN) ++fSumw2.fArray[bin];
|
---|
91 | if (biny == 0 || biny > fYaxis.GetNbins()) {
|
---|
92 | if (!fgStatOverflows) return -1;
|
---|
93 | }
|
---|
94 | ++fTsumw;
|
---|
95 | ++fTsumw2;
|
---|
96 | fTsumwy += y;
|
---|
97 | fTsumwy2 += y*y;
|
---|
98 | return bin;
|
---|
99 | }
|
---|
100 | ClassDef(MH2F, 1);
|
---|
101 | };
|
---|
102 |
|
---|
103 | class MProfile2D : public TProfile2D
|
---|
104 | {
|
---|
105 | public:
|
---|
106 | Int_t Fill(Double_t x, Double_t y, Double_t z)
|
---|
107 | {
|
---|
108 | Int_t bin,binx,biny;
|
---|
109 | fEntries++;
|
---|
110 | binx =TMath::Nint(x+1);
|
---|
111 | biny =fYaxis.FindBin(y);
|
---|
112 | bin = GetBin(binx, biny);
|
---|
113 | fArray[bin] += z;
|
---|
114 | fSumw2.fArray[bin] += z*z;
|
---|
115 | fBinEntries.fArray[bin] += 1;
|
---|
116 | if (fBinSumw2.fN) fBinSumw2.fArray[bin] += 1;
|
---|
117 | if (biny == 0 || biny > fYaxis.GetNbins()) {
|
---|
118 | if (!fgStatOverflows) return -1;
|
---|
119 | }
|
---|
120 | ++fTsumw;
|
---|
121 | ++fTsumw2;
|
---|
122 | fTsumwy += y;
|
---|
123 | fTsumwy2 += y*y;
|
---|
124 | fTsumwz += z;
|
---|
125 | fTsumwz2 += z*z;
|
---|
126 | return bin;
|
---|
127 | }
|
---|
128 | ClassDef(MProfile2D, 1);
|
---|
129 | };
|
---|
130 |
|
---|
131 |
|
---|
132 |
|
---|
133 | class MHSingles : public MH
|
---|
134 | {
|
---|
135 | MH2F fSignal;
|
---|
136 | MH2F fTime;
|
---|
137 | MProfile2D fPulse;
|
---|
138 |
|
---|
139 | UInt_t fNumEntries;
|
---|
140 |
|
---|
141 | MSingles *fSingles; //!
|
---|
142 | MPedestalSubtractedEvt *fEvent; //!
|
---|
143 |
|
---|
144 | public:
|
---|
145 | MHSingles() : fNumEntries(0), fSingles(0), fEvent(0)
|
---|
146 | {
|
---|
147 | fName = "MHSingles";
|
---|
148 |
|
---|
149 | fSignal.SetName("Signal");
|
---|
150 | fSignal.SetTitle("Title");
|
---|
151 | fSignal.SetXTitle("X title");
|
---|
152 | fSignal.SetYTitle("Y title");
|
---|
153 | fSignal.SetDirectory(NULL);
|
---|
154 |
|
---|
155 | fTime.SetName("Time");
|
---|
156 | fTime.SetTitle("Title");
|
---|
157 | fTime.SetXTitle("X title");
|
---|
158 | fTime.SetYTitle("Y title");
|
---|
159 | fTime.SetDirectory(NULL);
|
---|
160 |
|
---|
161 | fPulse.SetName("Pulse");
|
---|
162 | fPulse.SetTitle("Title");
|
---|
163 | fPulse.SetXTitle("X title");
|
---|
164 | fPulse.SetYTitle("Y title");
|
---|
165 | fPulse.SetDirectory(NULL);
|
---|
166 | }
|
---|
167 |
|
---|
168 | Bool_t SetupFill(const MParList *plist)
|
---|
169 | {
|
---|
170 | fSingles = (MSingles*)plist->FindObject("MSingles");
|
---|
171 | if (!fSingles)
|
---|
172 | {
|
---|
173 | *fLog << /* err << */ "MSingles not found... abort." << endl;
|
---|
174 | return kFALSE;
|
---|
175 | }
|
---|
176 |
|
---|
177 | fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
|
---|
178 | if (!fEvent)
|
---|
179 | {
|
---|
180 | *fLog << /* err << */ "MPedestalSubtractedEvt not found... abort." << endl;
|
---|
181 | return kFALSE;
|
---|
182 | }
|
---|
183 |
|
---|
184 | fNumEntries = 0;
|
---|
185 |
|
---|
186 | return kTRUE;
|
---|
187 | }
|
---|
188 | Bool_t ReInit(MParList *plist)
|
---|
189 | {
|
---|
190 | const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader");
|
---|
191 | if (!header)
|
---|
192 | {
|
---|
193 | *fLog << /* err << */ "MRawRunHeader not found... abort." << endl;
|
---|
194 | return kFALSE;
|
---|
195 | }
|
---|
196 |
|
---|
197 | const Int_t w = fSingles->GetIntegrationWindow();
|
---|
198 |
|
---|
199 | MBinning binsx, binsy;
|
---|
200 | binsx.SetEdges(fSingles->GetNumPixels(), -0.5, fSingles->GetNumPixels()-0.5);
|
---|
201 | binsy.SetEdges(300, -7.5*w, (60-7.5)*w);
|
---|
202 |
|
---|
203 | MH::SetBinning(fSignal, binsx, binsy);
|
---|
204 |
|
---|
205 | const UShort_t roi = header->GetNumSamples();
|
---|
206 |
|
---|
207 | // Correct for DRS timing!!
|
---|
208 | MBinning binst(roi, -0.5, roi-0.5);
|
---|
209 | MH::SetBinning(fTime, binsx, binst);
|
---|
210 |
|
---|
211 | MBinning binspy(2*roi, -roi-0.5, roi-0.5);
|
---|
212 | MH::SetBinning(fPulse, binsx, binspy);
|
---|
213 |
|
---|
214 | return kTRUE;
|
---|
215 | }
|
---|
216 |
|
---|
217 | Int_t Fill(const MParContainer *par, const Stat_t weight=1)
|
---|
218 | {
|
---|
219 | const Float_t *ptr = fEvent->GetSamples(0);
|
---|
220 | const Short_t roi = fEvent->GetNumSamples();
|
---|
221 |
|
---|
222 | for (unsigned int i=0; i<fSingles->GetNumPixels(); i++)
|
---|
223 | {
|
---|
224 | const vector<Single> &result = fSingles->GetVector(i);
|
---|
225 |
|
---|
226 | for (unsigned int j=0; j<result.size(); j++)
|
---|
227 | {
|
---|
228 | fSignal.Fill(i, result[j].fSignal);
|
---|
229 | fTime.Fill (i, result[j].fTime);
|
---|
230 |
|
---|
231 | if (!ptr)
|
---|
232 | continue;
|
---|
233 |
|
---|
234 | const Short_t tm = floor(result[j].fTime);
|
---|
235 |
|
---|
236 | for (int s=0; s<roi; s++)
|
---|
237 | fPulse.Fill(i, s-tm, ptr[s]);
|
---|
238 | }
|
---|
239 |
|
---|
240 | ptr += roi;
|
---|
241 | }
|
---|
242 |
|
---|
243 | fNumEntries++;
|
---|
244 |
|
---|
245 | return kTRUE;
|
---|
246 | }
|
---|
247 |
|
---|
248 | TH2 *GetSignal() { return &fSignal; }
|
---|
249 | TH2 *GetTime() { return &fTime; }
|
---|
250 | TH2 *GetPulse() { return &fPulse; }
|
---|
251 |
|
---|
252 | UInt_t GetNumEntries() const { return fNumEntries; }
|
---|
253 |
|
---|
254 | void Draw(Option_t *opt)
|
---|
255 | {
|
---|
256 | TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
|
---|
257 |
|
---|
258 | AppendPad("");
|
---|
259 |
|
---|
260 | pad->Divide(2,2);
|
---|
261 |
|
---|
262 | pad->cd(1);
|
---|
263 | gPad->SetBorderMode(0);
|
---|
264 | gPad->SetFrameBorderMode(0);
|
---|
265 | fSignal.Draw("colz");
|
---|
266 |
|
---|
267 | pad->cd(2);
|
---|
268 | gPad->SetBorderMode(0);
|
---|
269 | gPad->SetFrameBorderMode(0);
|
---|
270 | fTime.Draw("colz");
|
---|
271 |
|
---|
272 | pad->cd(3);
|
---|
273 | gPad->SetBorderMode(0);
|
---|
274 | gPad->SetFrameBorderMode(0);
|
---|
275 | fPulse.Draw("colz");
|
---|
276 | }
|
---|
277 |
|
---|
278 | ClassDef(MHSingles, 1)
|
---|
279 | };
|
---|
280 |
|
---|
281 | MStatusDisplay *d = new MStatusDisplay;
|
---|
282 |
|
---|
283 | MSingles *fSingles;
|
---|
284 |
|
---|
285 | TH1F *fluct1 = new TH1F("", "", 100, -50, 50);
|
---|
286 | TH1F *fluct2 = new TH1F("", "", 100, -50, 50);
|
---|
287 |
|
---|
288 | TGraph weights("weights.txt");
|
---|
289 |
|
---|
290 | Int_t PreProcess(MParList *plist)
|
---|
291 | {
|
---|
292 | fSingles = (MSingles*)plist->FindCreateObj("MSingles");
|
---|
293 | if (!fSingles)
|
---|
294 | return kFALSE;
|
---|
295 |
|
---|
296 | fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
|
---|
297 | if (!fEvent)
|
---|
298 | {
|
---|
299 | *fLog << /* err << */ "MPedestalSubtractedEvt not found... abort." << endl;
|
---|
300 | return kFALSE;
|
---|
301 | }
|
---|
302 |
|
---|
303 | d->AddTab("Fluct");
|
---|
304 | fluct2->Draw();
|
---|
305 | fluct1->Draw("same");
|
---|
306 |
|
---|
307 | fSingles->SetIntegrationWindow(20);
|
---|
308 |
|
---|
309 | //if (weights.GetN()==0)
|
---|
310 | // return kFALSE;
|
---|
311 |
|
---|
312 | return kTRUE;
|
---|
313 | }
|
---|
314 |
|
---|
315 | Int_t Process()
|
---|
316 | {
|
---|
317 | const UInt_t roi = fEvent->GetNumSamples();
|
---|
318 |
|
---|
319 | const size_t navg = 10;
|
---|
320 | const float threshold = 5;
|
---|
321 | const UInt_t start_skip = 20;
|
---|
322 | const UInt_t end_skip = 10;
|
---|
323 | const UInt_t integration_size = 10*3;
|
---|
324 |
|
---|
325 | vector<float> val(roi-navg);//result of first sliding average
|
---|
326 | for (int pix=0; pix<fEvent->GetNumPixels(); pix++)
|
---|
327 | {
|
---|
328 | vector<Single> &result = fSingles->GetVector(pix);
|
---|
329 | result.clear();
|
---|
330 |
|
---|
331 | const Float_t *ptr = fEvent->GetSamples(pix);
|
---|
332 |
|
---|
333 | //Sliding window
|
---|
334 | for (size_t i=0; i<roi-navg; i++)
|
---|
335 | {
|
---|
336 | val[i] = 0;
|
---|
337 | for (size_t j=i; j<i+navg; j++)
|
---|
338 | val[i] += ptr[j];
|
---|
339 | val[i] /= navg;
|
---|
340 | }
|
---|
341 |
|
---|
342 | //peak finding via threshold
|
---|
343 | for (UInt_t i=start_skip; i<roi-navg-end_skip-10; i++)
|
---|
344 | {
|
---|
345 | //search for threshold crossings
|
---|
346 | if (val[i+0]>threshold ||
|
---|
347 | val[i+4]>threshold ||
|
---|
348 |
|
---|
349 | val[i+5]<threshold ||
|
---|
350 | val[i+9]<threshold)
|
---|
351 | continue;
|
---|
352 |
|
---|
353 | // Evaluate arrival time more precisely!!!
|
---|
354 | // Get a better integration window
|
---|
355 |
|
---|
356 | Single single;
|
---|
357 | single.fSignal = 0;
|
---|
358 | single.fTime = i;
|
---|
359 |
|
---|
360 | // Crossing of the threshold is at 5
|
---|
361 | for (UInt_t j=i+5; j<i+5+integration_size; j++)
|
---|
362 | single.fSignal += ptr[j];
|
---|
363 |
|
---|
364 | result.push_back(single);
|
---|
365 |
|
---|
366 | i += 5+integration_size;
|
---|
367 | }
|
---|
368 | }
|
---|
369 | return kTRUE;
|
---|
370 | }
|
---|
371 |
|
---|
372 | Int_t PostProcess()
|
---|
373 | {
|
---|
374 | return kTRUE;
|
---|
375 | }
|
---|
376 |
|
---|
377 | Double_t fcn(Double_t *x, Double_t *par);
|
---|
378 |
|
---|
379 | int singlepe_final2()
|
---|
380 | {
|
---|
381 | // ======================================================
|
---|
382 |
|
---|
383 | const char *drsfile = "/fact/raw/2012/06/01/20120601_003.drs.fits.gz";
|
---|
384 |
|
---|
385 | MDirIter iter;
|
---|
386 | iter.AddDirectory("/fact/raw/2012/06/01", "20120601_005.fits.gz");
|
---|
387 |
|
---|
388 | // ======================================================
|
---|
389 |
|
---|
390 | // true: Display correctly mapped pixels in the camera displays
|
---|
391 | // but the value-vs-index plot is in software/spiral indices
|
---|
392 | // false: Display pixels in hardware/linear indices,
|
---|
393 | // but the order is the camera display is distorted
|
---|
394 | bool usemap = true;
|
---|
395 |
|
---|
396 | // map file to use (get that from La Palma!)
|
---|
397 | const char *map = usemap ? "/scratch_nfs/FACTmap111030.txt" : NULL;
|
---|
398 |
|
---|
399 | // ------------------------------------------------------
|
---|
400 |
|
---|
401 | Long_t max1 = 0;
|
---|
402 |
|
---|
403 | // ======================================================
|
---|
404 |
|
---|
405 | if (map && gSystem->AccessPathName(map, kFileExists))
|
---|
406 | {
|
---|
407 | gLog << "ERROR - Cannot access mapping file '" << map << "'" << endl;
|
---|
408 | return 1;
|
---|
409 | }
|
---|
410 | if (gSystem->AccessPathName(drsfile, kFileExists))
|
---|
411 | {
|
---|
412 | gLog << "ERROR - Cannot access drsfile file '" << drsfile << "'" << endl;
|
---|
413 | return 2;
|
---|
414 | }
|
---|
415 |
|
---|
416 | // --------------------------------------------------------------------------------
|
---|
417 |
|
---|
418 | gLog.Separator("Singles");
|
---|
419 | gLog << "Extract singles with '" << drsfile << "'" << endl;
|
---|
420 | gLog << endl;
|
---|
421 |
|
---|
422 | // ------------------------------------------------------
|
---|
423 |
|
---|
424 | MDrsCalibration drscalib300;
|
---|
425 | if (!drscalib300.ReadFits(drsfile))
|
---|
426 | return 4;
|
---|
427 |
|
---|
428 | // ------------------------------------------------------
|
---|
429 |
|
---|
430 | iter.Sort();
|
---|
431 | iter.Print();
|
---|
432 |
|
---|
433 | // ======================================================
|
---|
434 |
|
---|
435 | //MStatusDisplay *d = new MStatusDisplay;
|
---|
436 |
|
---|
437 | MBadPixelsCam badpixels;
|
---|
438 | badpixels.InitSize(1440);
|
---|
439 | badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
---|
440 | badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
---|
441 | badpixels[ 830].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
---|
442 | badpixels[ 923].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
---|
443 | badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
---|
444 | badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
---|
445 |
|
---|
446 | // Twin pixel
|
---|
447 | // 113
|
---|
448 | // 115
|
---|
449 | // 354
|
---|
450 | // 423
|
---|
451 | // 1195
|
---|
452 | // 1393
|
---|
453 |
|
---|
454 | //MDrsCalibrationTime timecam;
|
---|
455 |
|
---|
456 | MH::SetPalette();
|
---|
457 |
|
---|
458 | // ======================================================
|
---|
459 |
|
---|
460 | gLog << endl;
|
---|
461 | gLog.Separator("Processing external light pulser run");
|
---|
462 |
|
---|
463 | MTaskList tlist1;
|
---|
464 |
|
---|
465 | MSingles singles;
|
---|
466 | MRawRunHeader header;
|
---|
467 |
|
---|
468 | MParList plist1;
|
---|
469 | plist1.AddToList(&tlist1);
|
---|
470 | plist1.AddToList(&drscalib300);
|
---|
471 | plist1.AddToList(&badpixels);
|
---|
472 | plist1.AddToList(&singles);
|
---|
473 | plist1.AddToList(&header);
|
---|
474 |
|
---|
475 | MEvtLoop loop1("DetermineCalConst");
|
---|
476 | loop1.SetDisplay(d);
|
---|
477 | loop1.SetParList(&plist1);
|
---|
478 |
|
---|
479 | // ------------------ Setup the tasks ---------------
|
---|
480 |
|
---|
481 | MRawFitsRead read1;
|
---|
482 | read1.LoadMap(map);
|
---|
483 | read1.AddFiles(iter);
|
---|
484 |
|
---|
485 | MContinue cont1("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectCal");
|
---|
486 |
|
---|
487 | MGeomApply apply1;
|
---|
488 |
|
---|
489 | MDrsCalibApply drsapply1;
|
---|
490 |
|
---|
491 | MTaskInteractive mytask;
|
---|
492 | mytask.SetPreProcess(PreProcess);
|
---|
493 | mytask.SetProcess(Process);
|
---|
494 | mytask.SetPostProcess(PostProcess);
|
---|
495 |
|
---|
496 | MFillH fill("MHSingles");
|
---|
497 |
|
---|
498 | // ------------------ Setup eventloop and run analysis ---------------
|
---|
499 |
|
---|
500 | tlist1.AddToList(&read1);
|
---|
501 | // tlist1.AddToList(&cont1);
|
---|
502 | tlist1.AddToList(&apply1);
|
---|
503 | tlist1.AddToList(&drsapply1);
|
---|
504 | tlist1.AddToList(&mytask);
|
---|
505 | tlist1.AddToList(&fill);
|
---|
506 |
|
---|
507 | if (!loop1.Eventloop(max1))
|
---|
508 | return 10;
|
---|
509 |
|
---|
510 | if (!loop1.GetDisplay())
|
---|
511 | return 11;
|
---|
512 |
|
---|
513 | gStyle->SetOptFit(kTRUE);
|
---|
514 |
|
---|
515 | MHSingles* h = (MHSingles*) plist1.FindObject("MHSingles");
|
---|
516 | if (!h)
|
---|
517 | return 0;
|
---|
518 |
|
---|
519 | TH2 *hsignal = h->GetSignal();
|
---|
520 | TH2 *htime = h->GetTime();
|
---|
521 | TH2 *hpulse = h->GetPulse();
|
---|
522 |
|
---|
523 | d->AddTab("Time");
|
---|
524 | TH1 *ht = htime->ProjectionY();
|
---|
525 | ht->Scale(1./1440);
|
---|
526 | ht->Draw();
|
---|
527 |
|
---|
528 | d->AddTab("Pulse");
|
---|
529 | TH1 *hp = hpulse->ProjectionY();
|
---|
530 | hp->Scale(1./1440);
|
---|
531 | hp->Draw();
|
---|
532 |
|
---|
533 |
|
---|
534 | // ------------------ Spectrum Fit ---------------
|
---|
535 | // AFTER this the Code for the spektrum fit follows
|
---|
536 | // access the spektrumhistogram by the pointer *hist
|
---|
537 |
|
---|
538 | MGeomCamFACT fact;
|
---|
539 | MHCamera hcam(fact);
|
---|
540 |
|
---|
541 | //built an array of Amplitude spectra
|
---|
542 | TH1F hAmplitude("Rate", "Average number of dark counts per event",
|
---|
543 | 100, 0, 10);
|
---|
544 | TH1F hGain ("Gain", "Gain distribution",
|
---|
545 | 50, 100, 250);
|
---|
546 | TH1F hGainRMS ("RelSigma", "Distribution of rel. Sigma",
|
---|
547 | 100, 0, 30);
|
---|
548 | TH1F hCrosstlk ("Crosstalk", "Distribution of crosstalk probability",
|
---|
549 | 100, 0, 25);
|
---|
550 | TH1F hOffset ("Offset", "Baseline offset distribution",
|
---|
551 | 50, -7.5, 2.5);
|
---|
552 | TH1F hFitProb ("FitProb", "FitProb distribution",
|
---|
553 | 50, 0, 100);
|
---|
554 | TH1F hSum1 ("Sum1", "Sum",
|
---|
555 | 100, 0, 2000);
|
---|
556 | TH1F hSum2 ("Sum2", "Sum",
|
---|
557 | 100, 0, 10);
|
---|
558 |
|
---|
559 | // define fit function for Amplitudespectrum
|
---|
560 | TF1 func("spektrum", fcn, 0, 1000, 5);
|
---|
561 | func.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
|
---|
562 | func.SetLineColor(kBlue);
|
---|
563 |
|
---|
564 | // Map for which pixel shall be plotted and which not
|
---|
565 | TArrayC usePixel(1440);
|
---|
566 | memset(usePixel.GetArray(), 1, 1440);
|
---|
567 |
|
---|
568 | // List of Pixel that should be ignored in camera view
|
---|
569 | usePixel[424] = 0;
|
---|
570 | usePixel[923] = 0;
|
---|
571 | usePixel[1208] = 0;
|
---|
572 | usePixel[583] = 0;
|
---|
573 | usePixel[830] = 0;
|
---|
574 | usePixel[1399] = 0;
|
---|
575 | usePixel[113] = 0;
|
---|
576 | usePixel[115] = 0;
|
---|
577 | usePixel[354] = 0;
|
---|
578 | usePixel[423] = 0;
|
---|
579 | usePixel[1195] = 0;
|
---|
580 | usePixel[1393] = 0;
|
---|
581 |
|
---|
582 | int count_ok = 0;
|
---|
583 |
|
---|
584 | // Begin of Loop over Pixels
|
---|
585 | for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++)
|
---|
586 | {
|
---|
587 | if (usePixel[pixel]==0)
|
---|
588 | continue;
|
---|
589 |
|
---|
590 | //Projectipon of a certain Pixel out of Ampl.Spectrum
|
---|
591 | TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
|
---|
592 | hist->SetDirectory(0);
|
---|
593 |
|
---|
594 | const Int_t maxbin = hist->GetMaximumBin();
|
---|
595 | const Double_t binwidth = hist->GetBinWidth(maxbin);
|
---|
596 | const Double_t ampl = hist->GetBinContent(maxbin);
|
---|
597 |
|
---|
598 | double fwhm = 0;
|
---|
599 | for (int i=1; i<maxbin; i++)
|
---|
600 | if (hist->GetBinContent(maxbin-i)+hist->GetBinContent(maxbin+i)<ampl*2)
|
---|
601 | {
|
---|
602 | fwhm = (2*i+1)*binwidth;
|
---|
603 | break;
|
---|
604 | }
|
---|
605 |
|
---|
606 | if (fwhm==0)
|
---|
607 | {
|
---|
608 | cout << "Could not determine start value for sigma." << endl;
|
---|
609 | continue;
|
---|
610 | }
|
---|
611 |
|
---|
612 | // Amplitude, Gain, Sigma, Crosstalk Prob, Shift
|
---|
613 | // par[1] + par[4] is the position of the first maximum
|
---|
614 | Double_t par[5] =
|
---|
615 | {
|
---|
616 | ampl, hist->GetBinCenter(maxbin), fwhm*1.4, 0.1, -10
|
---|
617 | };
|
---|
618 |
|
---|
619 | func.SetParameters(par);
|
---|
620 |
|
---|
621 | const double min = par[1]-fwhm*1.4;
|
---|
622 | const double max = par[1]*5;
|
---|
623 |
|
---|
624 | //Fit Pixels spectrum
|
---|
625 | const TFitResultPtr rc = hist->Fit(&func, "N0QS", "", min, max);
|
---|
626 |
|
---|
627 | const bool ok = int(rc)==0;
|
---|
628 |
|
---|
629 | const double fit_prob = rc->Prob();
|
---|
630 | const float fGain = func.GetParameter(1);
|
---|
631 | const float fAmplitude = func.Integral(0, 7*fGain)/h->GetNumEntries(); //func.GetParameter(0);
|
---|
632 | const float f2GainRMS = func.GetParameter(2);
|
---|
633 | const float fCrosstlk = func.GetParameter(3);
|
---|
634 | const float fOffset = func.GetParameter(4)/singles.GetIntegrationWindow();
|
---|
635 |
|
---|
636 | hAmplitude.Fill(fAmplitude);
|
---|
637 | hGain.Fill(fGain);
|
---|
638 | //hGainRMS.Fill(f2GainRMS);
|
---|
639 | hCrosstlk.Fill(fCrosstlk*100);
|
---|
640 | hOffset.Fill(fOffset);
|
---|
641 | hFitProb.Fill(100*fit_prob);
|
---|
642 | hGainRMS.Fill(100*f2GainRMS/fGain);
|
---|
643 |
|
---|
644 | hcam.SetBinContent(pixel+1, fGain);
|
---|
645 |
|
---|
646 | usePixel[pixel] = ok;
|
---|
647 |
|
---|
648 | if (!ok)
|
---|
649 | {
|
---|
650 | cout << "Pixel " << setw(4) << pixel << ": failed!" << endl;
|
---|
651 | continue;
|
---|
652 | }
|
---|
653 |
|
---|
654 | for (int b=1; b<=hist->GetNbinsX(); b++)
|
---|
655 | {
|
---|
656 | hSum1.Fill(hist->GetBinCenter(b), hist->GetBinContent(b));
|
---|
657 | hSum2.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b));
|
---|
658 | }
|
---|
659 |
|
---|
660 | if (count_ok==0)
|
---|
661 | {
|
---|
662 | TCanvas &c = d->AddTab(Form("Pix%d", pixel));
|
---|
663 | c.cd();
|
---|
664 | hist->SetName(Form("Pix%d", pixel));
|
---|
665 | hist->DrawCopy()->SetDirectory(0);
|
---|
666 | //func.DrawCopy("same");
|
---|
667 | }
|
---|
668 |
|
---|
669 | count_ok++;
|
---|
670 |
|
---|
671 | delete hist;
|
---|
672 | }
|
---|
673 |
|
---|
674 | hcam.SetUsed(usePixel);
|
---|
675 |
|
---|
676 | TCanvas &canv = d->AddTab("Gain");
|
---|
677 | canv.cd();
|
---|
678 | hcam.SetNameTitle( "gain","Gain distribution over whole camera");
|
---|
679 | hcam.DrawCopy();
|
---|
680 |
|
---|
681 | TCanvas &c11 = d->AddTab("SumHist");
|
---|
682 | gPad->SetLogy();
|
---|
683 | hSum1.DrawCopy();
|
---|
684 |
|
---|
685 | TCanvas &c12 = d->AddTab("GainHist");
|
---|
686 | gPad->SetLogy();
|
---|
687 | hSum2.DrawCopy();
|
---|
688 |
|
---|
689 | TCanvas &c2 = d->AddTab("Result");
|
---|
690 | c2.Divide(3,2);
|
---|
691 | c2.cd(1);
|
---|
692 | gPad->SetLogy();
|
---|
693 | hAmplitude.DrawCopy();
|
---|
694 |
|
---|
695 | c2.cd(2);
|
---|
696 | gPad->SetLogy();
|
---|
697 | hGain.DrawCopy();
|
---|
698 |
|
---|
699 | c2.cd(3);
|
---|
700 | gPad->SetLogy();
|
---|
701 | hGainRMS.DrawCopy();
|
---|
702 |
|
---|
703 | c2.cd(4);
|
---|
704 | gPad->SetLogy();
|
---|
705 | hCrosstlk.DrawCopy();
|
---|
706 |
|
---|
707 | c2.cd(5);
|
---|
708 | gPad->SetLogy();
|
---|
709 | hOffset.DrawCopy();
|
---|
710 |
|
---|
711 | c2.cd(6);
|
---|
712 | gPad->SetLogy();
|
---|
713 | hFitProb.DrawCopy();
|
---|
714 |
|
---|
715 | /*
|
---|
716 | TCanvas &c3 = d->AddTab("Separation");
|
---|
717 | gPad->SetLogy();
|
---|
718 | hSep.DrawCopy();
|
---|
719 | */
|
---|
720 | return 0;
|
---|
721 | }
|
---|
722 |
|
---|
723 | Double_t fcn(Double_t *xx, Double_t *par)
|
---|
724 | {
|
---|
725 | Double_t x = xx[0];
|
---|
726 |
|
---|
727 | Double_t ampl = par[0];
|
---|
728 | Double_t gain = par[1];
|
---|
729 | Double_t sigma = par[2];
|
---|
730 | Double_t cross = par[3];
|
---|
731 | Double_t shift = par[4];
|
---|
732 |
|
---|
733 | Double_t xTalk = 1;
|
---|
734 |
|
---|
735 | Double_t y = 0;
|
---|
736 | for (int order = 1; order < 5; order++)
|
---|
737 | {
|
---|
738 | Double_t arg = (x - order*gain - shift)/sigma;
|
---|
739 |
|
---|
740 | Double_t gauss = TMath::Exp(-0.5 * arg * arg/order);
|
---|
741 |
|
---|
742 | y += xTalk*gauss;
|
---|
743 |
|
---|
744 | xTalk *= cross;
|
---|
745 | }
|
---|
746 |
|
---|
747 | return ampl*y;
|
---|
748 | }
|
---|