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 | #include <Getline.h>
|
---|
11 |
|
---|
12 | #include <cstdio>
|
---|
13 | #include <stdio.h>
|
---|
14 | #include <stdint.h>
|
---|
15 |
|
---|
16 | #include "MH.h"
|
---|
17 | #include "MLog.h"
|
---|
18 | #include "MLogManip.h"
|
---|
19 | #include "MDirIter.h"
|
---|
20 | #include "MFillH.h"
|
---|
21 | #include "MEvtLoop.h"
|
---|
22 | #include "MCamEvent.h"
|
---|
23 | #include "MGeomApply.h"
|
---|
24 | #include "MTaskList.h"
|
---|
25 | #include "MParList.h"
|
---|
26 | #include "MContinue.h"
|
---|
27 | #include "MBinning.h"
|
---|
28 | #include "MHDrsCalib.h"
|
---|
29 | #include "MDrsCalibApply.h"
|
---|
30 | #include "MRawFitsRead.h"
|
---|
31 | #include "MBadPixelsCam.h"
|
---|
32 | #include "MStatusDisplay.h"
|
---|
33 | #include "MTaskInteractive.h"
|
---|
34 | #include "MPedestalSubtractedEvt.h"
|
---|
35 | #include "MHCamera.h"
|
---|
36 | #include "MGeomCamFACT.h"
|
---|
37 | #include "MRawRunHeader.h"
|
---|
38 |
|
---|
39 | using namespace std;
|
---|
40 |
|
---|
41 | MPedestalSubtractedEvt *fEvent = 0;
|
---|
42 | MLog *fLog = &gLog;
|
---|
43 |
|
---|
44 | struct Single
|
---|
45 | {
|
---|
46 | float fSignal;
|
---|
47 | float fTime;
|
---|
48 | };
|
---|
49 |
|
---|
50 | class MSingles : public MParContainer, public MCamEvent
|
---|
51 | {
|
---|
52 | Int_t fIntegrationWindow;
|
---|
53 |
|
---|
54 | vector<vector<Single>> fData;
|
---|
55 |
|
---|
56 | public:
|
---|
57 | MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(0)
|
---|
58 | {
|
---|
59 | fName = name ? name : "MSingles";
|
---|
60 | }
|
---|
61 |
|
---|
62 | void InitSize(const UInt_t i)
|
---|
63 | {
|
---|
64 | fData.resize(i);
|
---|
65 | }
|
---|
66 |
|
---|
67 | vector<Single> &operator[](UInt_t i)
|
---|
68 | {
|
---|
69 | return fData[i];
|
---|
70 | }
|
---|
71 | vector<Single> &GetVector(UInt_t i)
|
---|
72 | {
|
---|
73 | return fData[i];
|
---|
74 | }
|
---|
75 |
|
---|
76 | UInt_t GetNumPixels() const
|
---|
77 | {
|
---|
78 | return fData.size();
|
---|
79 | }
|
---|
80 |
|
---|
81 | void SetIntegrationWindow(Int_t w)
|
---|
82 | {
|
---|
83 | fIntegrationWindow = w;
|
---|
84 | }
|
---|
85 | Int_t GetIntegrationWindow() const
|
---|
86 | {
|
---|
87 | return fIntegrationWindow;
|
---|
88 | }
|
---|
89 |
|
---|
90 | Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const
|
---|
91 | {
|
---|
92 | return kTRUE;
|
---|
93 | }
|
---|
94 | void DrawPixelContent(Int_t num) const { }
|
---|
95 |
|
---|
96 | ClassDef(MSingles, 1)
|
---|
97 | };
|
---|
98 |
|
---|
99 | class MH2F : public TH2F
|
---|
100 | {
|
---|
101 | public:
|
---|
102 | Int_t Fill(Double_t x,Double_t y)
|
---|
103 | {
|
---|
104 | Int_t binx, biny, bin;
|
---|
105 | fEntries++;
|
---|
106 | binx = TMath::Nint(x+1);
|
---|
107 | biny = fYaxis.FindBin(y);
|
---|
108 | bin = biny*(fXaxis.GetNbins()+2) + binx;
|
---|
109 | AddBinContent(bin);
|
---|
110 |
|
---|
111 | if (fSumw2.fN) ++fSumw2.fArray[bin];
|
---|
112 |
|
---|
113 | if (biny == 0 || biny > fYaxis.GetNbins())
|
---|
114 | {
|
---|
115 | if (!fgStatOverflows) return -1;
|
---|
116 | }
|
---|
117 |
|
---|
118 | ++fTsumw;
|
---|
119 | ++fTsumw2;
|
---|
120 | fTsumwy += y;
|
---|
121 | fTsumwy2 += y*y;
|
---|
122 | return bin;
|
---|
123 | }
|
---|
124 | ClassDef(MH2F, 1);
|
---|
125 | };
|
---|
126 |
|
---|
127 | class MProfile2D : public TProfile2D
|
---|
128 | {
|
---|
129 | public:
|
---|
130 | Int_t Fill(Double_t x, Double_t y, Double_t z)
|
---|
131 | {
|
---|
132 | Int_t bin,binx,biny;
|
---|
133 | fEntries++;
|
---|
134 | binx =TMath::Nint(x+1);
|
---|
135 | biny =fYaxis.FindBin(y);
|
---|
136 | bin = GetBin(binx, biny);
|
---|
137 | fArray[bin] += z;
|
---|
138 | fSumw2.fArray[bin] += z*z;
|
---|
139 | fBinEntries.fArray[bin] += 1;
|
---|
140 |
|
---|
141 | if (fBinSumw2.fN) fBinSumw2.fArray[bin] += 1;
|
---|
142 |
|
---|
143 | if (biny == 0 || biny > fYaxis.GetNbins())
|
---|
144 | {
|
---|
145 | if (!fgStatOverflows) return -1;
|
---|
146 | }
|
---|
147 |
|
---|
148 | ++fTsumw;
|
---|
149 | ++fTsumw2;
|
---|
150 | fTsumwy += y;
|
---|
151 | fTsumwy2 += y*y;
|
---|
152 | fTsumwz += z;
|
---|
153 | fTsumwz2 += z*z;
|
---|
154 | return bin;
|
---|
155 | }
|
---|
156 | ClassDef(MProfile2D, 1);
|
---|
157 | };
|
---|
158 |
|
---|
159 |
|
---|
160 |
|
---|
161 | class MHSingles : public MH
|
---|
162 | {
|
---|
163 | MH2F fSignal;
|
---|
164 | MH2F fTime;
|
---|
165 | MProfile2D fPulse;
|
---|
166 |
|
---|
167 | UInt_t fNumEntries;
|
---|
168 |
|
---|
169 | MSingles *fSingles; //!
|
---|
170 | MPedestalSubtractedEvt *fEvent; //!
|
---|
171 |
|
---|
172 | public:
|
---|
173 | MHSingles() : fNumEntries(0), fSingles(0), fEvent(0)
|
---|
174 | {
|
---|
175 | fName = "MHSingles";
|
---|
176 |
|
---|
177 | fSignal.SetName("Signal");
|
---|
178 | fSignal.SetTitle("pulse integral spectrum");
|
---|
179 | fSignal.SetXTitle("pixel [SoftID]");
|
---|
180 | fSignal.SetYTitle("time [a.u.]");
|
---|
181 | fSignal.SetDirectory(NULL);
|
---|
182 |
|
---|
183 | fTime.SetName("Time");
|
---|
184 | fTime.SetTitle("arival time spectrum");
|
---|
185 | fTime.SetXTitle("pixel [SoftID]");
|
---|
186 | fTime.SetYTitle("time [a.u.]");
|
---|
187 | fTime.SetDirectory(NULL);
|
---|
188 |
|
---|
189 | fPulse.SetName("Pulse");
|
---|
190 | fPulse.SetTitle("average pulse shape spectrum");
|
---|
191 | fPulse.SetXTitle("pixel [SoftID]");
|
---|
192 | fPulse.SetYTitle("time [a.u.]");
|
---|
193 | fPulse.SetDirectory(NULL);
|
---|
194 | }
|
---|
195 |
|
---|
196 | Bool_t SetupFill(const MParList *plist)
|
---|
197 | {
|
---|
198 | fSingles = (MSingles*)plist->FindObject("MSingles");
|
---|
199 |
|
---|
200 | if (!fSingles)
|
---|
201 | {
|
---|
202 | *fLog << /* err << */ "MSingles not found... abort." << endl;
|
---|
203 | return kFALSE;
|
---|
204 | }
|
---|
205 |
|
---|
206 | fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
|
---|
207 |
|
---|
208 | if (!fEvent)
|
---|
209 | {
|
---|
210 | *fLog << /* err << */ "MPedestalSubtractedEvt not found... abort." << endl;
|
---|
211 | return kFALSE;
|
---|
212 | }
|
---|
213 |
|
---|
214 | fNumEntries = 0;
|
---|
215 |
|
---|
216 | return kTRUE;
|
---|
217 | }
|
---|
218 | Bool_t ReInit(MParList *plist)
|
---|
219 | {
|
---|
220 | const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader");
|
---|
221 |
|
---|
222 | if (!header)
|
---|
223 | {
|
---|
224 | *fLog << /* err << */ "MRawRunHeader not found... abort." << endl;
|
---|
225 | return kFALSE;
|
---|
226 | }
|
---|
227 |
|
---|
228 | const Int_t w = fSingles->GetIntegrationWindow();
|
---|
229 |
|
---|
230 | MBinning binsx, binsy;
|
---|
231 |
|
---|
232 | binsx.SetEdges(fSingles->GetNumPixels(), -0.5, fSingles->GetNumPixels()-0.5);
|
---|
233 |
|
---|
234 | binsy.SetEdges(2*150, -2*7.5*w, 2*(60-7.5)*w);
|
---|
235 |
|
---|
236 | MH::SetBinning(fSignal, binsx, binsy);
|
---|
237 |
|
---|
238 | const UShort_t roi = header->GetNumSamples();
|
---|
239 |
|
---|
240 | // Correct for DRS timing!!
|
---|
241 | MBinning binst(roi, -0.5, roi-0.5);
|
---|
242 |
|
---|
243 | MH::SetBinning(fTime, binsx, binst);
|
---|
244 |
|
---|
245 | MBinning binspy(2*roi, -roi-0.5, roi-0.5);
|
---|
246 |
|
---|
247 | MH::SetBinning(fPulse, binsx, binspy);
|
---|
248 |
|
---|
249 | return kTRUE;
|
---|
250 | }
|
---|
251 |
|
---|
252 | Int_t Fill(const MParContainer *par, const Stat_t weight=1)
|
---|
253 | {
|
---|
254 | const Float_t *ptr = fEvent->GetSamples(0);
|
---|
255 | const Short_t roi = fEvent->GetNumSamples();
|
---|
256 |
|
---|
257 | for (unsigned int i=0; i<fSingles->GetNumPixels(); i++)
|
---|
258 | {
|
---|
259 | const vector<Single> &result = fSingles->GetVector(i);
|
---|
260 |
|
---|
261 | for (unsigned int j=0; j<result.size(); j++)
|
---|
262 | {
|
---|
263 | fSignal.Fill(i, result[j].fSignal);
|
---|
264 | fTime.Fill (i, result[j].fTime);
|
---|
265 |
|
---|
266 | if (!ptr)
|
---|
267 | continue;
|
---|
268 |
|
---|
269 | const Short_t tm = floor(result[j].fTime);
|
---|
270 |
|
---|
271 | for (int s=0; s<roi; s++)
|
---|
272 | fPulse.Fill(i, s-tm, ptr[s]);
|
---|
273 | }
|
---|
274 |
|
---|
275 | ptr += roi;
|
---|
276 | }
|
---|
277 |
|
---|
278 | fNumEntries++;
|
---|
279 |
|
---|
280 | return kTRUE;
|
---|
281 | }
|
---|
282 |
|
---|
283 | TH2 *GetSignal()
|
---|
284 | {
|
---|
285 | return &fSignal;
|
---|
286 | }
|
---|
287 | TH2 *GetTime()
|
---|
288 | {
|
---|
289 | return &fTime;
|
---|
290 | }
|
---|
291 | TH2 *GetPulse()
|
---|
292 | {
|
---|
293 | return &fPulse;
|
---|
294 | }
|
---|
295 |
|
---|
296 | UInt_t GetNumEntries() const
|
---|
297 | {
|
---|
298 | return fNumEntries;
|
---|
299 | }
|
---|
300 |
|
---|
301 | void Draw(Option_t *opt)
|
---|
302 | {
|
---|
303 | TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
|
---|
304 |
|
---|
305 | AppendPad("");
|
---|
306 |
|
---|
307 | pad->Divide(2,2);
|
---|
308 |
|
---|
309 | pad->cd(1);
|
---|
310 | gPad->SetBorderMode(0);
|
---|
311 | gPad->SetFrameBorderMode(0);
|
---|
312 | fSignal.Draw("colz");
|
---|
313 |
|
---|
314 | pad->cd(2);
|
---|
315 | gPad->SetBorderMode(0);
|
---|
316 | gPad->SetFrameBorderMode(0);
|
---|
317 | fTime.Draw("colz");
|
---|
318 |
|
---|
319 | pad->cd(3);
|
---|
320 | gPad->SetBorderMode(0);
|
---|
321 | gPad->SetFrameBorderMode(0);
|
---|
322 | fPulse.Draw("colz");
|
---|
323 | }
|
---|
324 |
|
---|
325 | ClassDef(MHSingles, 1)
|
---|
326 | };
|
---|
327 |
|
---|
328 | MStatusDisplay *d = new MStatusDisplay;
|
---|
329 |
|
---|
330 | MSingles *fSingles;
|
---|
331 |
|
---|
332 | TH1F *fluct1 = new TH1F("", "", 200, -200, 200);
|
---|
333 | TH1F *fluct2 = new TH1F("", "", 100, -50, 50);
|
---|
334 |
|
---|
335 | TGraph weights("weights.txt");
|
---|
336 |
|
---|
337 | Int_t PreProcess(MParList *plist)
|
---|
338 | {
|
---|
339 | fSingles = (MSingles*)plist->FindCreateObj("MSingles");
|
---|
340 |
|
---|
341 | if (!fSingles)
|
---|
342 | return kFALSE;
|
---|
343 |
|
---|
344 | fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
|
---|
345 |
|
---|
346 | if (!fEvent)
|
---|
347 | {
|
---|
348 | *fLog << /* err << */ "MPedestalSubtractedEvt not found... abort." << endl;
|
---|
349 | return kFALSE;
|
---|
350 | }
|
---|
351 |
|
---|
352 | d->AddTab("Fluct");
|
---|
353 | fluct2->Draw();
|
---|
354 | fluct1->Draw("same");
|
---|
355 |
|
---|
356 | fSingles->SetIntegrationWindow(20);
|
---|
357 |
|
---|
358 | //if (weights.GetN()==0)
|
---|
359 | // return kFALSE;
|
---|
360 |
|
---|
361 | return kTRUE;
|
---|
362 | }
|
---|
363 |
|
---|
364 | Int_t Process()
|
---|
365 | {
|
---|
366 | const UInt_t roi = fEvent->GetNumSamples();
|
---|
367 |
|
---|
368 | const size_t navg = 10;
|
---|
369 | const float threshold = 5;
|
---|
370 | const UInt_t start_skip = 20;
|
---|
371 | const UInt_t end_skip = 10;
|
---|
372 | const UInt_t integration_size = 10*3;
|
---|
373 | const UInt_t max_search_window = 30;
|
---|
374 |
|
---|
375 | vector<float> val(roi-navg);//result of first sliding average
|
---|
376 |
|
---|
377 | for (int pix=0; pix<fEvent->GetNumPixels(); pix++)
|
---|
378 | {
|
---|
379 | vector<Single> &result = fSingles->GetVector(pix);
|
---|
380 | result.clear();
|
---|
381 |
|
---|
382 | const Float_t *ptr = fEvent->GetSamples(pix);
|
---|
383 |
|
---|
384 | //Sliding window
|
---|
385 | for (size_t i=0; i<roi-navg; i++)
|
---|
386 | {
|
---|
387 | val[i] = 0;
|
---|
388 |
|
---|
389 | for (size_t j=i; j<i+navg; j++)
|
---|
390 | val[i] += ptr[j];
|
---|
391 |
|
---|
392 | val[i] /= navg;
|
---|
393 | }
|
---|
394 |
|
---|
395 | //peak finding via threshold
|
---|
396 | for (UInt_t i=start_skip; i<roi-navg-end_skip-30; i++)
|
---|
397 | {
|
---|
398 | //search for threshold crossings
|
---|
399 | if (val[i+0]>threshold ||
|
---|
400 | val[i+4]>threshold ||
|
---|
401 |
|
---|
402 | val[i+5]<threshold ||
|
---|
403 | val[i+9]<threshold)
|
---|
404 | continue;
|
---|
405 |
|
---|
406 | //search for maximum after threshold crossing
|
---|
407 | UInt_t k_max = i+5;
|
---|
408 |
|
---|
409 | for (UInt_t k=i+5; k<i+max_search_window; k++)
|
---|
410 | {
|
---|
411 | if (val[k] > val[k_max])
|
---|
412 | k_max = k;
|
---|
413 | }
|
---|
414 |
|
---|
415 | if (k_max == i+5 || k_max == i + max_search_window-1) continue;
|
---|
416 |
|
---|
417 | //search for half maximum before maximum
|
---|
418 | UInt_t k_half_max = 0;
|
---|
419 |
|
---|
420 | for (UInt_t k=k_max; k>k_max-25; k--)
|
---|
421 | {
|
---|
422 | if (val[k-1] < val[k_max]/2 &&
|
---|
423 | val[k] >= val[k_max]/2 )
|
---|
424 | {
|
---|
425 | k_half_max = k;
|
---|
426 | break;
|
---|
427 | }
|
---|
428 | }
|
---|
429 |
|
---|
430 | if (k_half_max > roi-navg-end_skip-max_search_window - integration_size)
|
---|
431 | continue;
|
---|
432 |
|
---|
433 | if (k_half_max == 0)
|
---|
434 | continue;
|
---|
435 |
|
---|
436 | if (k_max - k_half_max > 14)
|
---|
437 | continue;
|
---|
438 |
|
---|
439 | fluct2->Fill(k_max - k_half_max);
|
---|
440 |
|
---|
441 | // Evaluate arrival time more precisely!!!
|
---|
442 | // Get a better integration window
|
---|
443 |
|
---|
444 | Single single;
|
---|
445 | single.fSignal = 0;
|
---|
446 | single.fTime = k_half_max;
|
---|
447 |
|
---|
448 |
|
---|
449 | // Crossing of the threshold is at 5
|
---|
450 | for (UInt_t j=k_half_max; j<k_half_max+integration_size; j++)
|
---|
451 | single.fSignal += ptr[j];
|
---|
452 |
|
---|
453 | result.push_back(single);
|
---|
454 |
|
---|
455 | i += 5+integration_size;
|
---|
456 | }
|
---|
457 | }
|
---|
458 |
|
---|
459 | return kTRUE;
|
---|
460 | }
|
---|
461 |
|
---|
462 | Int_t PostProcess()
|
---|
463 | {
|
---|
464 | return kTRUE;
|
---|
465 | }
|
---|
466 |
|
---|
467 | Double_t fcn(Double_t *x, Double_t *par);
|
---|
468 | Double_t fcn_cross(Double_t *x, Double_t *par);
|
---|
469 |
|
---|
470 | void FitGaussiansToSpectrum( UInt_t nfunc, TH1 &hSource, TH1 &hDest,
|
---|
471 | Int_t maxbin, double fwhm, Double_t gain );
|
---|
472 |
|
---|
473 | Double_t MedianOfH1 (
|
---|
474 | TH1* inputHisto
|
---|
475 | );
|
---|
476 |
|
---|
477 | int singlepe(
|
---|
478 | // const char *runfile, const char *drsfile, const char *outfile
|
---|
479 | )
|
---|
480 | {
|
---|
481 |
|
---|
482 | // ======================================================
|
---|
483 |
|
---|
484 | const char *drsfile = "/fact/raw/2012/05/18/20120518_003.drs.fits.gz";
|
---|
485 | // const char *drsfile = "/fact/raw/2012/06/25/20120625_112.drs.fits.gz";
|
---|
486 |
|
---|
487 | MDirIter iter;
|
---|
488 | iter.AddDirectory("/fact/raw/2012/05/18", "20120518_005.fits.gz");
|
---|
489 | // iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(runfile));
|
---|
490 |
|
---|
491 | // TString filename;
|
---|
492 | // for (UInt_t fileNr = 116; fileNr < 147; fileNr++)
|
---|
493 | // {
|
---|
494 | // filename = "20120625_";
|
---|
495 | // filename += fileNr;
|
---|
496 | // filename += ".fits.gz";
|
---|
497 | // iter.AddDirectory("/fact/raw/2012/06/25", filename);
|
---|
498 | // }
|
---|
499 |
|
---|
500 | // iter.AddDirectory("/fact/raw/2012/06/25", "20120625_116.fits.gz");
|
---|
501 | // iter.AddDirectory("/fact/raw/2012/06/25", "20120625_117.fits.gz");
|
---|
502 | // iter.AddDirectory("/fact/raw/2012/06/25", "20120625_118.fits.gz");
|
---|
503 | // iter.AddDirectory("/fact/raw/2012/06/25", "20120625_119.fits.gz");
|
---|
504 |
|
---|
505 | // ======================================================
|
---|
506 |
|
---|
507 | // true: Display correctly mapped pixels in the camera displays
|
---|
508 | // but the value-vs-index plot is in software/spiral indices
|
---|
509 | // false: Display pixels in hardware/linear indices,
|
---|
510 | // but the order is the camera display is distorted
|
---|
511 | bool usemap = true;
|
---|
512 |
|
---|
513 | // map file to use (get that from La Palma!)
|
---|
514 | const char *map = usemap ? "/scratch_nfs/FACTmap111030.txt" : NULL;
|
---|
515 |
|
---|
516 | // ------------------------------------------------------
|
---|
517 |
|
---|
518 | Long_t max1 = 0;
|
---|
519 |
|
---|
520 | // ======================================================
|
---|
521 |
|
---|
522 | if (map && gSystem->AccessPathName(map, kFileExists))
|
---|
523 | {
|
---|
524 | gLog << "ERROR - Cannot access mapping file '" << map << "'" << endl;
|
---|
525 | return 1;
|
---|
526 | }
|
---|
527 |
|
---|
528 | if (gSystem->AccessPathName(drsfile, kFileExists))
|
---|
529 | {
|
---|
530 | gLog << "ERROR - Cannot access drsfile file '" << drsfile << "'" << endl;
|
---|
531 | return 2;
|
---|
532 | }
|
---|
533 |
|
---|
534 | // --------------------------------------------------------------------------------
|
---|
535 |
|
---|
536 | gLog.Separator("Singles");
|
---|
537 | gLog << "Extract singles with '" << drsfile << "'" << endl;
|
---|
538 | gLog << endl;
|
---|
539 |
|
---|
540 | // ------------------------------------------------------
|
---|
541 |
|
---|
542 | MDrsCalibration drscalib300;
|
---|
543 |
|
---|
544 | if (!drscalib300.ReadFits(drsfile))
|
---|
545 | return 4;
|
---|
546 |
|
---|
547 | // ------------------------------------------------------
|
---|
548 |
|
---|
549 | iter.Sort();
|
---|
550 | iter.Print();
|
---|
551 |
|
---|
552 | // ======================================================
|
---|
553 |
|
---|
554 | //MStatusDisplay *d = new MStatusDisplay;
|
---|
555 |
|
---|
556 | MBadPixelsCam badpixels;
|
---|
557 | badpixels.InitSize(1440);
|
---|
558 | badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
---|
559 | badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
---|
560 | badpixels[ 830].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
---|
561 | badpixels[ 923].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
---|
562 | badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
---|
563 | badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable);
|
---|
564 |
|
---|
565 | // Twin pixel
|
---|
566 | // 113
|
---|
567 | // 115
|
---|
568 | // 354
|
---|
569 | // 423
|
---|
570 | // 1195
|
---|
571 | // 1393
|
---|
572 |
|
---|
573 | //MDrsCalibrationTime timecam;
|
---|
574 |
|
---|
575 | MH::SetPalette();
|
---|
576 |
|
---|
577 | // ======================================================
|
---|
578 |
|
---|
579 | gLog << endl;
|
---|
580 | gLog.Separator("Processing external light pulser run");
|
---|
581 |
|
---|
582 | MTaskList tlist1;
|
---|
583 |
|
---|
584 | MSingles singles;
|
---|
585 | MRawRunHeader header;
|
---|
586 |
|
---|
587 | MParList plist1;
|
---|
588 | plist1.AddToList(&tlist1);
|
---|
589 | plist1.AddToList(&drscalib300);
|
---|
590 | plist1.AddToList(&badpixels);
|
---|
591 | plist1.AddToList(&singles);
|
---|
592 | plist1.AddToList(&header);
|
---|
593 |
|
---|
594 | TString Title;
|
---|
595 | Title = iter.Next();
|
---|
596 | iter.Reset();
|
---|
597 | Title += "; ";
|
---|
598 | Title += max1;
|
---|
599 |
|
---|
600 | MEvtLoop loop1(Title);
|
---|
601 | loop1.SetDisplay(d);
|
---|
602 | loop1.SetParList(&plist1);
|
---|
603 |
|
---|
604 | // ------------------ Setup the tasks ---------------
|
---|
605 |
|
---|
606 | MRawFitsRead read1;
|
---|
607 | read1.LoadMap(map);
|
---|
608 | read1.AddFiles(iter);
|
---|
609 |
|
---|
610 | MContinue cont1("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectCal");
|
---|
611 |
|
---|
612 | MGeomApply apply1;
|
---|
613 |
|
---|
614 | MDrsCalibApply drsapply1;
|
---|
615 |
|
---|
616 | MTaskInteractive mytask;
|
---|
617 | mytask.SetPreProcess(PreProcess);
|
---|
618 | mytask.SetProcess(Process);
|
---|
619 | mytask.SetPostProcess(PostProcess);
|
---|
620 |
|
---|
621 | MFillH fill("MHSingles");
|
---|
622 |
|
---|
623 | // ------------------ Setup eventloop and run analysis ---------------
|
---|
624 |
|
---|
625 | tlist1.AddToList(&read1);
|
---|
626 | // tlist1.AddToList(&cont1);
|
---|
627 | tlist1.AddToList(&apply1);
|
---|
628 | tlist1.AddToList(&drsapply1);
|
---|
629 | tlist1.AddToList(&mytask);
|
---|
630 | tlist1.AddToList(&fill);
|
---|
631 |
|
---|
632 | if (!loop1.Eventloop(max1))
|
---|
633 | return 10;
|
---|
634 |
|
---|
635 | if (!loop1.GetDisplay())
|
---|
636 | return 11;
|
---|
637 |
|
---|
638 | gStyle->SetOptFit(kTRUE);
|
---|
639 |
|
---|
640 | MHSingles* h = (MHSingles*) plist1.FindObject("MHSingles");
|
---|
641 |
|
---|
642 | if (!h)
|
---|
643 | return 0;
|
---|
644 |
|
---|
645 | TH2 *hsignal = h->GetSignal();
|
---|
646 | TH2 *htime = h->GetTime();
|
---|
647 | TH2 *hpulse = h->GetPulse();
|
---|
648 |
|
---|
649 | d->AddTab("Time");
|
---|
650 | TH1 *ht = htime->ProjectionY();
|
---|
651 | ht->SetYTitle("Counts [a.u.]");
|
---|
652 | ht->Scale(1./1440);
|
---|
653 | ht->Draw();
|
---|
654 |
|
---|
655 | d->AddTab("Pulse");
|
---|
656 | TH1 *hp = hpulse->ProjectionY();
|
---|
657 | hp->SetYTitle("Counts [a.u.]");
|
---|
658 | hp->Scale(1./1440);
|
---|
659 | hp->Draw();
|
---|
660 |
|
---|
661 |
|
---|
662 | // ------------------ Spectrum Fit ---------------
|
---|
663 | // AFTER this the Code for the spektrum fit follows
|
---|
664 | // access the spektrumhistogram by the pointer *hist
|
---|
665 | UInt_t maxOrder = 8;
|
---|
666 |
|
---|
667 |
|
---|
668 | MGeomCamFACT fact;
|
---|
669 | MHCamera hcam(fact);
|
---|
670 | MHCamera hcam2(fact);
|
---|
671 |
|
---|
672 | //built an array of Amplitude spectra
|
---|
673 | TH1F hAmplitude ("Rate", "Average number of dark counts per event",
|
---|
674 | 200, 0, 20);
|
---|
675 | hAmplitude.SetXTitle("Amplitude [a.u.]");
|
---|
676 | hAmplitude.SetYTitle("Rate");
|
---|
677 |
|
---|
678 | TH1F hGain ("Gain", "Gain distribution",
|
---|
679 | 400, 100, 300);
|
---|
680 | hGain.SetXTitle("gain [a.u.]");
|
---|
681 | hGain.SetYTitle("Rate");
|
---|
682 |
|
---|
683 | TH1F hGainRMS ("RelSigma", "Distribution of rel. Sigma",
|
---|
684 | 100, 0, 30);
|
---|
685 | hGainRMS.SetXTitle("gainRMS [a.u.]");
|
---|
686 | hGainRMS.SetYTitle("Rate");
|
---|
687 | hGainRMS.SetStats(false);
|
---|
688 |
|
---|
689 | TH1F hCrosstlk ("Crosstalk", "Distribution of crosstalk probability",
|
---|
690 | 100, 0, 25);
|
---|
691 | hCrosstlk.SetXTitle("crosstalk [a.u.]");
|
---|
692 | hCrosstlk.SetYTitle("Rate");
|
---|
693 |
|
---|
694 | TH1F hOffset ("Offset", "Baseline offset distribution",
|
---|
695 | 50, -7.5, 2.5);
|
---|
696 | hOffset.SetXTitle("baseline offset [a.u.]");
|
---|
697 | hOffset.SetYTitle("Rate");
|
---|
698 |
|
---|
699 | TH1F hFitProb ("FitProb", "FitProb distribution",
|
---|
700 | 50, 0, 100);
|
---|
701 | hFitProb.SetXTitle("FitProb p-value [a.u.]");
|
---|
702 | hFitProb.SetYTitle("Rate");
|
---|
703 | hFitProb.SetStats(false);
|
---|
704 |
|
---|
705 |
|
---|
706 | TH1F hSum ("Sum1", "Sum of all pixel's' integral spectra",
|
---|
707 | 100, 0, 2000);
|
---|
708 | hSum.SetXTitle("pulse integral [a.u.]");
|
---|
709 | hSum.SetYTitle("Rate");
|
---|
710 | hSum.SetStats(false);
|
---|
711 | hSum.Sumw2();
|
---|
712 |
|
---|
713 |
|
---|
714 | TH1F hSumScale ("Sum2",
|
---|
715 | "Sum of all pixel's' integral spectra (scaled with gain)",
|
---|
716 | 100, 0.05, 10.05);
|
---|
717 | hSumScale.SetXTitle("pulse integral [pe]");
|
---|
718 | hSumScale.SetYTitle("Rate");
|
---|
719 | hSumScale.SetStats(false);
|
---|
720 | hSumScale.Sumw2();
|
---|
721 |
|
---|
722 | // define fit function for Amplitudespectrum
|
---|
723 | TF1 func("spektrum", fcn, 0, 1000, 5);
|
---|
724 | func.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
|
---|
725 | func.SetLineColor(kBlue);
|
---|
726 |
|
---|
727 | // define fit function for Amplitudespectrum
|
---|
728 | TF1 funcSum("sum_spektrum", fcn, 0, 2000, 5);
|
---|
729 | funcSum.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset", "PowOfXtalk");
|
---|
730 | funcSum.SetLineColor(kBlue);
|
---|
731 |
|
---|
732 | // define fit function for Amplitudespectrum
|
---|
733 | TF1 funcScaled("gain_sum_spektrum", fcn, 0, 10, 5);
|
---|
734 | funcScaled.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
|
---|
735 | funcScaled.SetLineColor(kBlue);
|
---|
736 |
|
---|
737 | TF1 funcScaledBL("gain_sum_spektrum", fcn, 0, 10, 5);
|
---|
738 | funcScaledBL.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
|
---|
739 | funcScaledBL.SetLineColor(kRed);
|
---|
740 |
|
---|
741 |
|
---|
742 |
|
---|
743 | // Map for which pixel shall be plotted and which not
|
---|
744 | TArrayC usePixel(1440);
|
---|
745 | memset(usePixel.GetArray(), 1, 1440);
|
---|
746 |
|
---|
747 | // List of Pixel that should be ignored in camera view
|
---|
748 | usePixel[424] = 0;
|
---|
749 | usePixel[923] = 0;
|
---|
750 | usePixel[1208] = 0;
|
---|
751 | usePixel[583] = 0;
|
---|
752 | usePixel[830] = 0;
|
---|
753 | usePixel[1399] = 0;
|
---|
754 | usePixel[113] = 0;
|
---|
755 | usePixel[115] = 0;
|
---|
756 | usePixel[354] = 0;
|
---|
757 | usePixel[423] = 0;
|
---|
758 | usePixel[1195] = 0;
|
---|
759 | usePixel[1393] = 0;
|
---|
760 |
|
---|
761 | // usePixel[990] = 0;
|
---|
762 |
|
---|
763 | // Map for which pixel which were suspicous
|
---|
764 | bool suspicous[1440];
|
---|
765 |
|
---|
766 | for (int i=0; i<1440; i++) suspicous[i]=false;
|
---|
767 |
|
---|
768 | // List of Pixel that should be ignored in camera view
|
---|
769 | // suspicous[802] = true;
|
---|
770 |
|
---|
771 | //------------------------------------------------------------------------
|
---|
772 | // Fill and calculate sum spectrum
|
---|
773 |
|
---|
774 | // Begin of Loop over Pixels
|
---|
775 | for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++)
|
---|
776 | {
|
---|
777 | //jump to next pixel if the current is marked as faulty
|
---|
778 | if (usePixel[pixel]==0)
|
---|
779 | continue;
|
---|
780 |
|
---|
781 | TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
|
---|
782 | hist->SetDirectory(0);
|
---|
783 |
|
---|
784 | //loop over slices
|
---|
785 | for (int b=1; b<=hist->GetNbinsX(); b++)
|
---|
786 | {
|
---|
787 | //Fill sum spectrum
|
---|
788 | hSum.Fill(hist->GetBinCenter(b), hist->GetBinContent(b));
|
---|
789 | }
|
---|
790 | }
|
---|
791 |
|
---|
792 | for (UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
|
---|
793 | {
|
---|
794 | hSum.SetBinError(bin, hSum.GetBinContent(bin)*0.1);
|
---|
795 | }
|
---|
796 |
|
---|
797 | gROOT->SetSelectedPad(0);
|
---|
798 | TCanvas &c11 = d->AddTab("SumHist");
|
---|
799 | c11.cd();
|
---|
800 | gPad->SetLogy();
|
---|
801 | gPad->SetGridx();
|
---|
802 | gPad->SetGridy();
|
---|
803 | hSum.DrawCopy("hist");
|
---|
804 | //------------------------fit sum spectrum-------------------------------
|
---|
805 | const Int_t maxbin = hSum.GetMaximumBin();
|
---|
806 | const Double_t binwidth = hSum.GetBinWidth(maxbin);
|
---|
807 | const Double_t ampl = hSum.GetBinContent(maxbin);
|
---|
808 |
|
---|
809 | double fwhmSum = 0;
|
---|
810 |
|
---|
811 | //Calculate full width half Maximum
|
---|
812 | for (int i=1; i<maxbin; i++)
|
---|
813 | if (hSum.GetBinContent(maxbin-i)+hSum.GetBinContent(maxbin+i)<ampl*2)
|
---|
814 | {
|
---|
815 | fwhmSum = (2*i+1)*binwidth;
|
---|
816 | break;
|
---|
817 | }
|
---|
818 |
|
---|
819 | if (fwhmSum==0)
|
---|
820 | {
|
---|
821 | cout << "Could not determine start value for sigma." << endl;
|
---|
822 | }
|
---|
823 |
|
---|
824 | //Set fit parameters
|
---|
825 | Double_t sum_par[6] =
|
---|
826 | {
|
---|
827 | ampl, hSum.GetBinCenter(maxbin), fwhmSum*1.4, 0.1, -10, 1
|
---|
828 | };
|
---|
829 | funcSum.SetParameters(sum_par);
|
---|
830 |
|
---|
831 | //calculate fit range
|
---|
832 | const double min = sum_par[1]-fwhmSum*2;
|
---|
833 | const double max = sum_par[1]*(maxOrder);
|
---|
834 | funcSum.SetRange(min, max);
|
---|
835 | funcSum.FixParameter(5,0.1);
|
---|
836 |
|
---|
837 | //Fit and draw spectrum
|
---|
838 | hSum.Fit(&funcSum, "NOQS", "", min, max);
|
---|
839 | funcSum.DrawCopy("SAME");
|
---|
840 | funcSum.GetParameters(sum_par);
|
---|
841 |
|
---|
842 | //Set Parameters for following fits
|
---|
843 | // const float Amplitude = sum_par[0];
|
---|
844 | const float gain = sum_par[1];
|
---|
845 | // const float GainRMS = sum_par[2];
|
---|
846 | const float Crosstlk = sum_par[3];
|
---|
847 | const float Offset = sum_par[4];
|
---|
848 | // const float PowCrosstlk = sum_par[5];
|
---|
849 |
|
---|
850 |
|
---|
851 | //--------fit gausses to peaks of sum specetrum -----------
|
---|
852 |
|
---|
853 | // create histo for height of Maximum amplitudes vs. pe
|
---|
854 | double min_bin = hSum.GetBinCenter(maxbin);
|
---|
855 | double max_bin = hSum.GetBinCenter((maxOrder-2)*(maxbin));
|
---|
856 | double bin_width= (max_bin - min_bin)/10;
|
---|
857 |
|
---|
858 | TH1D hMaxHeightsSum("MaxHeightSum", "peak heights",
|
---|
859 | 10, min_bin-bin_width, max_bin*2-bin_width/2 );
|
---|
860 |
|
---|
861 | FitGaussiansToSpectrum
|
---|
862 | (
|
---|
863 | maxOrder-2,
|
---|
864 | hSum,
|
---|
865 | hMaxHeightsSum,
|
---|
866 | maxbin,
|
---|
867 | fwhmSum,
|
---|
868 | gain
|
---|
869 | );
|
---|
870 |
|
---|
871 | //EOF: fit sum spectrum-----------------------------
|
---|
872 |
|
---|
873 | hMaxHeightsSum.SetLineColor(kRed);
|
---|
874 | // hMaxHeightsSum.DrawCopy("SAME");
|
---|
875 |
|
---|
876 | //Fit Peak heights
|
---|
877 | TF1 fExpo1( "fExpo1","expo", min, max);
|
---|
878 | fExpo1.SetLineColor(kRed);
|
---|
879 | hMaxHeightsSum.Fit(&fExpo1, "N0QS" );
|
---|
880 | // hMaxHeightsSum.DrawCopy("SAME");
|
---|
881 | // fExpo1.DrawCopy("SAME");
|
---|
882 |
|
---|
883 | // EOF: Fill and calculate sum spectrum
|
---|
884 | //------------------------------------------------------------------------
|
---|
885 |
|
---|
886 |
|
---|
887 |
|
---|
888 | //------------------------------------------------------------------------
|
---|
889 | // Gain Calculation for each Pixel
|
---|
890 |
|
---|
891 | int count_ok = 0;
|
---|
892 |
|
---|
893 | // Begin of Loop over Pixels
|
---|
894 | for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++)
|
---|
895 | {
|
---|
896 |
|
---|
897 | if (usePixel[pixel]==0)
|
---|
898 | continue;
|
---|
899 |
|
---|
900 | //Projectipon of a certain Pixel out of Ampl.Spectrum
|
---|
901 | TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
|
---|
902 | hist->SetDirectory(0);
|
---|
903 |
|
---|
904 | for (UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
|
---|
905 | {
|
---|
906 | hist->SetBinError(bin, hSum.GetBinContent(bin)*0.1);
|
---|
907 | }
|
---|
908 |
|
---|
909 | const Int_t maxbin = hist->GetMaximumBin();
|
---|
910 |
|
---|
911 | const Double_t binwidth = hist->GetBinWidth(maxbin);
|
---|
912 |
|
---|
913 | const Double_t ampl = hist->GetBinContent(maxbin);
|
---|
914 |
|
---|
915 | const Double_t maxPos = hist->GetBinCenter(maxbin);
|
---|
916 |
|
---|
917 | double fwhm = 0;
|
---|
918 |
|
---|
919 | for (int i=1; i<maxbin; i++)
|
---|
920 | if (hist->GetBinContent(maxbin-i)+hist->GetBinContent(maxbin+i)<ampl*2)
|
---|
921 | {
|
---|
922 | fwhm = (2*i+1)*binwidth;
|
---|
923 | break;
|
---|
924 | }
|
---|
925 |
|
---|
926 | if (fwhm==0)
|
---|
927 | {
|
---|
928 | cout << "Could not determine start value for sigma." << endl;
|
---|
929 | continue;
|
---|
930 | }
|
---|
931 |
|
---|
932 | // Amplitude, Gain, Sigma, Crosstalk Prob, Shift
|
---|
933 | // par[1] + par[4] is the position of the first maximum
|
---|
934 | Double_t par[5] =
|
---|
935 | {
|
---|
936 | ampl,
|
---|
937 | hist->GetBinCenter(maxbin),
|
---|
938 | // gain,
|
---|
939 | fwhm*2,
|
---|
940 | // GainRMS,
|
---|
941 | Crosstlk,
|
---|
942 | Offset
|
---|
943 | };
|
---|
944 |
|
---|
945 |
|
---|
946 | func.SetParameters(par);
|
---|
947 | const double fit_min = par[1]-fwhm*1.4;
|
---|
948 | const double fit_max = par[1]*maxOrder;
|
---|
949 | func.SetRange(fit_min-fwhm*2, fit_max);
|
---|
950 | func.SetParLimits(1, maxPos - fwhm, maxPos + fwhm);
|
---|
951 |
|
---|
952 | if (suspicous[pixel] == true)
|
---|
953 | {
|
---|
954 | cout << "Maxbin\t" << maxbin << endl;
|
---|
955 | cout << "min\t" << fit_min << endl;
|
---|
956 | cout << "max\t" << fit_max << endl;
|
---|
957 | cout << "Amplitude, 1.Max, Sigma, fwhm:\t"
|
---|
958 | << par[0] << "\t"
|
---|
959 | << par[1] << "\t"
|
---|
960 | << par[2] << "\t"
|
---|
961 | << fwhm << "\t" << endl;
|
---|
962 | }
|
---|
963 |
|
---|
964 | //Rebin Projection
|
---|
965 | hist->Rebin(2);
|
---|
966 |
|
---|
967 | //Fit Pixels spectrum
|
---|
968 | const TFitResultPtr rc = hist->Fit(&func, "N0QS", "", fit_min, fit_max);
|
---|
969 |
|
---|
970 | const bool ok = int(rc)==0;
|
---|
971 |
|
---|
972 | const double fit_prob = rc->Prob();
|
---|
973 | const float fGain = func.GetParameter(1);
|
---|
974 | const float fAmplitude = func.Integral(0, 7*fGain)/h->GetNumEntries(); //func.GetParameter(0);
|
---|
975 | const float f2GainRMS = func.GetParameter(2);
|
---|
976 | const float fCrosstlk = func.GetParameter(3);
|
---|
977 | const float fOffset = func.GetParameter(4)/singles.GetIntegrationWindow();
|
---|
978 | // const float fCrossOrd = func.GetParameter(5);
|
---|
979 |
|
---|
980 | hAmplitude.Fill(fAmplitude);
|
---|
981 | hGain.Fill(fGain);
|
---|
982 | //hGainRMS.Fill(f2GainRMS);
|
---|
983 | hCrosstlk.Fill(fCrosstlk*100);
|
---|
984 | hOffset.Fill(fOffset);
|
---|
985 | hFitProb.Fill(100*fit_prob);
|
---|
986 | hGainRMS.Fill(100*f2GainRMS/fGain);
|
---|
987 |
|
---|
988 | hcam.SetBinContent(pixel+1, fGain);
|
---|
989 | hcam2.SetBinContent(pixel+1, 100*f2GainRMS/fGain);
|
---|
990 |
|
---|
991 | usePixel[pixel] = ok;
|
---|
992 |
|
---|
993 | if (!ok)
|
---|
994 | {
|
---|
995 | cout << "Pixel " << setw(4) << pixel << ": failed!" << endl;
|
---|
996 | continue;
|
---|
997 | }
|
---|
998 |
|
---|
999 | //Fill sum spectrum
|
---|
1000 | for (int b=1; b<=hist->GetNbinsX(); b++)
|
---|
1001 | {
|
---|
1002 | // hSum.Fill(hist->GetBinCenter(b), hist->GetBinContent(b));
|
---|
1003 | hSumScale.Fill((hist->GetBinCenter(b)-fOffset)/fGain, (hist->GetBinContent(b)));
|
---|
1004 | }
|
---|
1005 |
|
---|
1006 | //plott special pixel
|
---|
1007 | if (
|
---|
1008 | count_ok==0 ||
|
---|
1009 | suspicous[pixel]
|
---|
1010 | )
|
---|
1011 | {
|
---|
1012 | TCanvas &c = d->AddTab(Form("Pix%d", pixel));
|
---|
1013 | c.cd();
|
---|
1014 | gPad->SetLogy();
|
---|
1015 | gPad->SetGridx();
|
---|
1016 | gPad->SetGridy();
|
---|
1017 |
|
---|
1018 | hist->SetStats(false);
|
---|
1019 | hist->SetYTitle("Counts [a.u.]");
|
---|
1020 | hist->SetName(Form("Pix%d", pixel));
|
---|
1021 | hist->DrawCopy("hist")->SetDirectory(0);
|
---|
1022 | // hist->Draw();
|
---|
1023 | func.DrawCopy("SAME");
|
---|
1024 | }
|
---|
1025 |
|
---|
1026 | count_ok++;
|
---|
1027 |
|
---|
1028 | delete hist;
|
---|
1029 | // if (suspicous[pixel] == true) usePixel[pixel]=0;
|
---|
1030 | }
|
---|
1031 |
|
---|
1032 | // End of Loop over Pixel
|
---|
1033 |
|
---|
1034 | //------------------fill histograms-----------------------
|
---|
1035 |
|
---|
1036 | hcam.SetUsed(usePixel);
|
---|
1037 | hcam2.SetUsed(usePixel);
|
---|
1038 |
|
---|
1039 | TCanvas &canv = d->AddTab("Gain");
|
---|
1040 | canv.cd();
|
---|
1041 | canv.Divide(2,2);
|
---|
1042 |
|
---|
1043 | canv.cd(1);
|
---|
1044 | hcam.SetNameTitle( "gain","Gain distribution over whole camera");
|
---|
1045 | hcam.DrawCopy();
|
---|
1046 |
|
---|
1047 | canv.cd(2);
|
---|
1048 | hcam2.SetNameTitle( "gainRMS","GainRMS distribution over whole camera");
|
---|
1049 | hcam2.DrawCopy();
|
---|
1050 |
|
---|
1051 | TCanvas &cam_canv = d->AddTab("Camera_Gain");
|
---|
1052 |
|
---|
1053 | //------------------ Draw gain corrected sum specetrum -------------------
|
---|
1054 | gROOT->SetSelectedPad(0);
|
---|
1055 | TCanvas &c12 = d->AddTab("GainHist");
|
---|
1056 | c12.cd();
|
---|
1057 | gPad->SetLogy();
|
---|
1058 | gPad->SetGridx();
|
---|
1059 | gPad->SetGridy();
|
---|
1060 |
|
---|
1061 | for (UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
|
---|
1062 | {
|
---|
1063 | hSumScale.SetBinError(bin, hSum.GetBinContent(bin)*0.1);
|
---|
1064 | }
|
---|
1065 |
|
---|
1066 | hSumScale.DrawCopy("hist");
|
---|
1067 |
|
---|
1068 |
|
---|
1069 | //-------------------- fit gain corrected sum spectrum -------------------
|
---|
1070 | const Int_t maxbin2 = hSumScale.GetMaximumBin();
|
---|
1071 | const Double_t binwidth2 = hSumScale.GetBinWidth(maxbin2);
|
---|
1072 | const Double_t ampl2 = hSumScale.GetBinContent(maxbin2);
|
---|
1073 |
|
---|
1074 | //Calculate full width half Maximum
|
---|
1075 | double fwhmScaled = 0;
|
---|
1076 |
|
---|
1077 | for (int i=1; i<maxbin; i++)
|
---|
1078 | if (hSumScale.GetBinContent(maxbin2-i)+hSumScale.GetBinContent(maxbin2+i)<ampl2*2)
|
---|
1079 | {
|
---|
1080 | fwhmScaled = (2*i+1)*binwidth2;
|
---|
1081 | break;
|
---|
1082 | }
|
---|
1083 |
|
---|
1084 | if (fwhmScaled==0)
|
---|
1085 | {
|
---|
1086 | cout << "Could not determine start value for sigma." << endl;
|
---|
1087 | }
|
---|
1088 |
|
---|
1089 | //Set fit parameters
|
---|
1090 | Double_t par2[6] =
|
---|
1091 | {
|
---|
1092 | ampl2, 0.9, fwhmScaled*1.4, Crosstlk, 0, 1.1
|
---|
1093 | };
|
---|
1094 |
|
---|
1095 | funcScaled.SetParameters(par2);
|
---|
1096 | funcScaledBL.SetParameters(par2);
|
---|
1097 | // funcScaled.FixParameter(1,0.9);
|
---|
1098 | funcScaled.FixParameter(4,0);
|
---|
1099 | funcScaledBL.FixParameter(4,0);
|
---|
1100 | funcScaled.FixParameter(5,Crosstlk);
|
---|
1101 | funcScaledBL.FixParameter(5,Crosstlk);
|
---|
1102 |
|
---|
1103 | const double minScaled = par2[1]-fwhmScaled;
|
---|
1104 | const double maxScaled = par2[1]*maxOrder;
|
---|
1105 | cout << "par2[1]" << par2[1] << endl;
|
---|
1106 | cout << "maxScaled\t" << maxScaled
|
---|
1107 | << " \t\t minScaled\t" << minScaled << endl;
|
---|
1108 |
|
---|
1109 | funcScaled.SetRange(minScaled, maxScaled);
|
---|
1110 | funcScaledBL.SetRange(minScaled, maxScaled);
|
---|
1111 | hSumScale.Fit(&funcScaled, "N0QS", "", minScaled, maxScaled);
|
---|
1112 | hSumScale.Fit(&funcScaledBL, "WLN0QS", "", minScaled, maxScaled);
|
---|
1113 | funcScaled.DrawCopy("SAME");
|
---|
1114 | funcScaledBL.DrawCopy("SAME");
|
---|
1115 | //--------fit gausses to peaks of gain corrected sum specetrum -----------
|
---|
1116 |
|
---|
1117 | // create histo for height of Maximum amplitudes vs. pe
|
---|
1118 | TH1D hMaxHeights("MaxHeights", "peak heights",
|
---|
1119 | 10, hSumScale.GetBinCenter(maxbin-1)-0.5, 10+hSumScale.GetBinCenter(maxbin-1)-0.5);
|
---|
1120 |
|
---|
1121 | FitGaussiansToSpectrum
|
---|
1122 | (
|
---|
1123 | maxOrder-2,
|
---|
1124 | hSumScale,
|
---|
1125 | hMaxHeights,
|
---|
1126 | maxbin2,
|
---|
1127 | fwhmScaled,
|
---|
1128 | 1
|
---|
1129 | );
|
---|
1130 |
|
---|
1131 |
|
---|
1132 |
|
---|
1133 | // TCanvas &c16 = d->AddTab("PeakHeights");
|
---|
1134 | // gPad->SetLogy();
|
---|
1135 | hMaxHeights.SetLineColor(kRed);
|
---|
1136 | // hMaxHeights.DrawCopy("SAME");
|
---|
1137 | TF1 fExpo( "fExpo","expo", 0, maxOrder-2);
|
---|
1138 | fExpo.SetLineColor(kRed);
|
---|
1139 | hMaxHeights.Fit(&fExpo, "N0QS" );
|
---|
1140 | // fExpo.DrawCopy("SAME");
|
---|
1141 | // c12.cd();
|
---|
1142 | // fExpo.DrawCopy("SAME");
|
---|
1143 |
|
---|
1144 | TCanvas &c2 = d->AddTab("Result");
|
---|
1145 | c2.Divide(3,2);
|
---|
1146 | c2.cd(1);
|
---|
1147 | gPad->SetLogy();
|
---|
1148 | hAmplitude.DrawCopy();
|
---|
1149 |
|
---|
1150 | c2.cd(2);
|
---|
1151 | gPad->SetLogy();
|
---|
1152 | // hGain.DrawCopy();
|
---|
1153 |
|
---|
1154 |
|
---|
1155 | TF1 GainGaus( "GainGaus", "gaus", 0, hGain.GetXaxis()->GetXmax());
|
---|
1156 | hGain.Rebin(2);
|
---|
1157 | hGain.Fit(&GainGaus, "N0QS");
|
---|
1158 |
|
---|
1159 | // float gain_mean = GainGaus.GetParameter(1);
|
---|
1160 | float gain_median = MedianOfH1(&hGain);
|
---|
1161 | TH1F hNormGain ("NormGain", "median normalized Gain distribution",
|
---|
1162 | hGain.GetNbinsX(),
|
---|
1163 | (hGain.GetXaxis()->GetXmin())/gain_median,
|
---|
1164 | (hGain.GetXaxis()->GetXmax())/gain_median
|
---|
1165 | );
|
---|
1166 | hNormGain.SetXTitle("gain [fraction of median gain value]");
|
---|
1167 | hNormGain.SetYTitle("Counts");
|
---|
1168 |
|
---|
1169 |
|
---|
1170 | hNormGain.Add(&hGain);
|
---|
1171 | // hNormGain.SetStats(false);
|
---|
1172 | hNormGain.GetXaxis()->SetRangeUser(
|
---|
1173 | 1-(hNormGain.GetXaxis()->GetXmax()-1),
|
---|
1174 | hNormGain.GetXaxis()->GetXmax()
|
---|
1175 | );
|
---|
1176 | hNormGain.DrawCopy();
|
---|
1177 | hNormGain.Fit(&GainGaus, "N0QS");
|
---|
1178 | GainGaus.DrawCopy("SAME");
|
---|
1179 |
|
---|
1180 | //plot normailzed camera view
|
---|
1181 | cam_canv.cd();
|
---|
1182 | MHCamera hNormCam(fact);
|
---|
1183 | hNormCam.SetStats(false);
|
---|
1184 |
|
---|
1185 | for (UInt_t pixel =1; pixel <= 1440; pixel++)
|
---|
1186 | {
|
---|
1187 | hNormCam.SetBinContent(pixel, hcam.GetBinContent(pixel)/gain_median);
|
---|
1188 | }
|
---|
1189 |
|
---|
1190 | hNormCam.SetNameTitle( "gain","normalized [Median] Gain distribution over whole camera");
|
---|
1191 | hNormCam.SetZTitle("Horst");
|
---|
1192 | hNormCam.SetUsed(usePixel);
|
---|
1193 | hNormCam.DrawCopy();
|
---|
1194 | // hGain.Scale(1/GainGaus.GetParameter(1));
|
---|
1195 | // GainGaus.DrawCopy("SAME");
|
---|
1196 |
|
---|
1197 | c2.cd(3);
|
---|
1198 | gPad->SetLogy();
|
---|
1199 | hGainRMS.DrawCopy();
|
---|
1200 |
|
---|
1201 | c2.cd(4);
|
---|
1202 | gPad->SetLogy();
|
---|
1203 | hCrosstlk.DrawCopy();
|
---|
1204 |
|
---|
1205 | c2.cd(5);
|
---|
1206 | gPad->SetLogy();
|
---|
1207 | hOffset.DrawCopy();
|
---|
1208 |
|
---|
1209 | c2.cd(6);
|
---|
1210 | gPad->SetLogy();
|
---|
1211 | // hFitProb.DrawCopy();
|
---|
1212 | hGain.DrawCopy();
|
---|
1213 |
|
---|
1214 |
|
---|
1215 | TCanvas &c25 = d->AddTab("Spectra");
|
---|
1216 | c25.Divide(2,1);
|
---|
1217 | c25.cd(1);
|
---|
1218 | gPad->SetLogy();
|
---|
1219 | gPad->SetGrid();
|
---|
1220 | hSum.DrawCopy("hist");
|
---|
1221 | funcSum.DrawCopy("SAME");
|
---|
1222 |
|
---|
1223 | c25.cd(2);
|
---|
1224 | gPad->SetLogy();
|
---|
1225 | gPad->SetGrid();
|
---|
1226 | hSumScale.DrawCopy("hist");
|
---|
1227 | funcScaled.DrawCopy("SAME");
|
---|
1228 | funcScaledBL.DrawCopy("SAME");
|
---|
1229 |
|
---|
1230 |
|
---|
1231 | /*
|
---|
1232 | TCanvas &c3 = d->AddTab("Separation");
|
---|
1233 | gPad->SetLogy();
|
---|
1234 | hSep.DrawCopy();
|
---|
1235 | */
|
---|
1236 |
|
---|
1237 | d->SaveAs("/home_nfs/isdc/jbbuss/singlepe.root");
|
---|
1238 | return 0;
|
---|
1239 | }
|
---|
1240 |
|
---|
1241 | Double_t fcn(Double_t *xx, Double_t *par)
|
---|
1242 | {
|
---|
1243 | Double_t x = xx[0];
|
---|
1244 |
|
---|
1245 | Double_t ampl = par[0];
|
---|
1246 | Double_t gain = par[1];
|
---|
1247 | Double_t sigma = par[2];
|
---|
1248 | Double_t cross = par[3];
|
---|
1249 | Double_t shift = par[4];
|
---|
1250 | // Double_t k = par[5];
|
---|
1251 |
|
---|
1252 | Double_t xTalk = 1;
|
---|
1253 |
|
---|
1254 | Double_t y = 0;
|
---|
1255 |
|
---|
1256 | for (int order = 1; order < 14; order++)
|
---|
1257 | {
|
---|
1258 | Double_t arg = (x - order*gain - shift)/sigma;
|
---|
1259 |
|
---|
1260 | Double_t gauss = TMath::Exp(-0.5 * arg * arg/order);
|
---|
1261 |
|
---|
1262 | y += xTalk*gauss;
|
---|
1263 |
|
---|
1264 | // for (int j = 1; j < k; j++)
|
---|
1265 | // {
|
---|
1266 | xTalk *= cross;
|
---|
1267 | // }
|
---|
1268 | }
|
---|
1269 |
|
---|
1270 | return ampl*y;
|
---|
1271 | }
|
---|
1272 |
|
---|
1273 | Double_t fcn_cross(Double_t *xx, Double_t *par)
|
---|
1274 | {
|
---|
1275 | Double_t x = xx[0];
|
---|
1276 |
|
---|
1277 | Double_t ampl = par[0];
|
---|
1278 | Double_t gain = par[1];
|
---|
1279 | Double_t sigma = par[2];
|
---|
1280 | Double_t cross = par[3];
|
---|
1281 | Double_t shift = par[4];
|
---|
1282 | // Double_t powOfX = par[5];
|
---|
1283 |
|
---|
1284 | Double_t xTalk = 1;
|
---|
1285 |
|
---|
1286 | Double_t y = 0;
|
---|
1287 |
|
---|
1288 | for (int order = 1; order < 14; order++)
|
---|
1289 | {
|
---|
1290 | Double_t arg = (x - order*gain - shift)/sigma;
|
---|
1291 |
|
---|
1292 | Double_t gauss = TMath::Exp(-0.5 * arg * arg/order);
|
---|
1293 |
|
---|
1294 | y += xTalk*gauss;
|
---|
1295 |
|
---|
1296 | // xTalk = pow(cross, order*powOfX);
|
---|
1297 | // for (int j = 1; j < powOfX; j++)
|
---|
1298 | // {
|
---|
1299 | // xTalk *= cross;
|
---|
1300 | // }
|
---|
1301 | // cross *= TMath::Exp(-powOfX*cross);
|
---|
1302 | xTalk *= cross;
|
---|
1303 | }
|
---|
1304 |
|
---|
1305 | return ampl*y;
|
---|
1306 | }
|
---|
1307 |
|
---|
1308 |
|
---|
1309 | void
|
---|
1310 | FitGaussiansToSpectrum(
|
---|
1311 | UInt_t maxOrder,
|
---|
1312 | TH1 &hSource,
|
---|
1313 | TH1 &hDest,
|
---|
1314 | Int_t maxbin,
|
---|
1315 | double fwhm,
|
---|
1316 | Double_t gain
|
---|
1317 | )
|
---|
1318 | {
|
---|
1319 |
|
---|
1320 | Double_t peakposition = hSource.GetBinCenter(maxbin);
|
---|
1321 | char fname[20];
|
---|
1322 |
|
---|
1323 | // fit gauss functions to spectrum
|
---|
1324 | for (UInt_t nr = 1; nr<maxOrder; nr++)
|
---|
1325 | {
|
---|
1326 | sprintf(fname,"gaus%d",nr);
|
---|
1327 |
|
---|
1328 | //create gauss functions
|
---|
1329 | TF1 gaus( fname,"gaus", peakposition-fwhm, peakposition+fwhm);
|
---|
1330 | gaus.SetParLimits(1, peakposition-fwhm, peakposition+fwhm);
|
---|
1331 | gaus.SetParLimits(2, -fwhm, fwhm );
|
---|
1332 | //fit and draw gauss functions
|
---|
1333 | hSource.Fit( &gaus, "N0QS", "", peakposition-fwhm, peakposition+fwhm);
|
---|
1334 | // gaus.DrawCopy("SAME");
|
---|
1335 |
|
---|
1336 | peakposition = gaus.GetParameter(1);
|
---|
1337 |
|
---|
1338 | //fill histo for height of Maximum amplitudes vs. pe
|
---|
1339 | hDest.SetBinContent(nr, gaus.GetParameter(0));
|
---|
1340 |
|
---|
1341 | peakposition += gain;
|
---|
1342 | }
|
---|
1343 |
|
---|
1344 | return;
|
---|
1345 | }
|
---|
1346 |
|
---|
1347 | Double_t MedianOfH1 (
|
---|
1348 | TH1* inputHisto
|
---|
1349 | )
|
---|
1350 | {
|
---|
1351 | //compute the median for 1-d histogram h1
|
---|
1352 | Int_t nbins = inputHisto->GetXaxis()->GetNbins();
|
---|
1353 | Double_t *x = new Double_t[nbins];
|
---|
1354 | Double_t *y = new Double_t[nbins];
|
---|
1355 |
|
---|
1356 | for (Int_t i=0; i<nbins; i++)
|
---|
1357 | {
|
---|
1358 | x[i] = inputHisto->GetXaxis()->GetBinCenter(i+1);
|
---|
1359 | y[i] = inputHisto->GetBinContent(i+1);
|
---|
1360 | }
|
---|
1361 |
|
---|
1362 | Double_t median = TMath::Median(nbins,x,y);
|
---|
1363 | delete [] x;
|
---|
1364 | delete [] y;
|
---|
1365 | return median;
|
---|
1366 | }
|
---|
1367 | // end of PlotMedianEachSliceOfPulse
|
---|
1368 | //----------------------------------------------------------------------------
|
---|