bug fix
[u/mrichter/AliRoot.git] / FMD / scripts / FindCommonModeNoise.C
CommitLineData
506dc39d 1//____________________________________________________________________
2//
3// $Id: DrawDigits.C 30718 2009-01-22 16:07:40Z cholm $
4//
5// Script that contains a class to draw eloss from hits, versus ADC
6// counts from digits, using the AliFMDInputHits class in the util library.
7//
8// It draws the energy loss versus the p/(mq^2). It can be overlayed
9// with the Bethe-Bloc curve to show how the simulation behaves
10// relative to the expected.
11//
12// Use the script `Compile.C' to compile this class using ACLic.
13//
14#include <TH1F.h>
15#include <TMath.h>
16#include <TStyle.h>
17#include <TArrayF.h>
18#include <TSystem.h>
19#include <TFile.h>
20#include <TCanvas.h>
21#include <TSystem.h>
22#include <TRandom.h>
23#include <TPad.h>
24#include <TLegend.h>
25
26#include <AliCDBManager.h>
27#include <AliFMDDigit.h>
28#include <AliFMDInput.h>
29#include <AliFMDParameters.h>
30#include <AliRawReader.h>
31#include <AliLog.h>
32
33#include <iostream>
34
35/** @class DrawDigits
36 @brief Draw pedestal, noise, and both common mode noise corrected
37 @code
38 Root> .L Compile.C
39 Root> Compile("FindCommonModeNoise.C")
40 Root> FindCommonModeNoise fcmn("file");
41 Root> fcmn.Run();
42 @endcode
43
44 If null is passed instead of a file, then a simulation of common
45 mode noise is run instead.
46
47 @ingroup FMD_script
48 */
49class FindCommonModeNoise : public AliFMDInput
50{
51private:
52 TString fCalibDir;
53 Int_t fCount;
54 TFile* fOut;
55 Float_t fShift;
56 Float_t fCmn; // Common mode noise component
57 TH1F* fPed; //
58 TH1F* fNoise; //
59 TH1F* fCenteredPed; //
60 TH1F* fCorrectedNoise; //
61 TH1F* fSummedAdc; //
62 TH1F* fSummedAdc2; //
63 TH1F* fSummedCenteredAdc; //
64 TH1F* fSummedCenteredAdc2; //
65 TH1F* fAdc; //
66
67public:
68 //__________________________________________________________________
69 FindCommonModeNoise(const char* file,
70 const char* calibDir=0)
71 : AliFMDInput(),
72 fCalibDir(""),
73 fOut(0)
74 {
75 fCmn = 2;
76 if (calibDir) fCalibDir = calibDir;
77 if (file) {
78 AddLoad(kRaw);
79 SetRawFile(file);
80 }
81 else
82 AddLoad(kUser);
83
84 TString outName(Form("histo_%s", gSystem->BaseName(fRawFile.Data())));
85 if (!outName.EndsWith(".root")) outName.Append(".root");
86 fOut = TFile::Open(outName.Data(),"RECREATE");
87
88 fPed = MakeHist("ped_dist","Pedestal Distribution");
89 fNoise = MakeHist("noise_dist",
90 "Noise Distribution (Not CMN corrected)");
91 fCenteredPed = MakeHist("cmn_ped_dist",
92 "Pedestal - Average Pedestal Distribution");
93 fCorrectedNoise = MakeHist("noise_corrected_dist",
94 "Noise Distribution (CMN corrected)");
95 fSummedAdc = MakeHist("sum_adc_dist","Sum ADC Distribution");
96 fSummedAdc2 = MakeHist("sum_adc2_dist",
97 "Sum ADC Squared Distribution");
98 fSummedCenteredAdc = MakeHist("sum_cmn_adc_dist",
99 "Sum CMN Adjusted ADC Distribution");
100 fSummedCenteredAdc2 = MakeHist("sum_cmn_adc2_dist",
101 "Sum CMN Adjusted ADC Squared Distribution");
102 fAdc = MakeHist("adc_dist","Current ADC Distribution");
103 }
104 //__________________________________________________________________
105 Int_t NEvents() const
106 {
107 Int_t nEv = (fReader ? fReader->GetNumberOfEvents() : -1);
108 if (nEv > 0) return TMath::Min(nEv, 1000);
109 return 1000;
110 }
111 //__________________________________________________________________
112 TH1F* MakeHist(const char* name, const char* title)
113 {
114 TH1F* h = new TH1F(name,title,51200,-0.5,51199.5);
115 h->SetXTitle("Strip number");
116 h->SetFillColor(kRed);
117 h->SetFillStyle(3001);
118 h->SetDirectory(0);
119 h->SetStats(kFALSE);
120 return h;
121 }
122 //__________________________________________________________________
123 Bool_t Init()
124 {
125 AliCDBManager* cdb = AliCDBManager::Instance();
126 cdb->SetRun(0);
127 cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
128
129 AliFMDParameters* pars = AliFMDParameters::Instance();
130 if (fCalibDir.IsNull())
131 pars->Init(kFALSE);
132 else
133 pars->Init(fCalibDir.Data(), kFALSE);
134
135 fCount = 0;
136
137 return AliFMDInput::Init();
138 }
139
140 //__________________________________________________________________
141 Bool_t Begin(Int_t e)
142 {
143 fCount++;
144 fAdc->Reset();
145 // Common mode noise generation
146 fShift = gRandom->Gaus(0, fCmn);
147 return AliFMDInput::Begin(e);
148 }
149
150 //__________________________________________________________________
151 Int_t BinCenter(UShort_t d, Char_t r, UShort_t s, UShort_t t) const
152 {
153 if (d == 1 && !(r == 'I' || r == 'i')) return -1;
154
155 Int_t idx = -1;
156 switch (d) {
157 case 1: idx = 0; break;
158 case 2: idx = 10240; break;
159 case 3: idx = 3 * 10240; break;
160 default:
161 return -1;
162 }
163
164 UShort_t q = (r == 'I' || r == 'i' ? 0 : 1);
165 // UShort_t nSec = (q == 0 ? 20 : 40);
166 UShort_t nStr = (q == 0 ? 512 : 256);
167
168 idx += (q == 0 ? 0 : 10240);
169 idx += s * nStr + t;
170
171 if (idx >= 51200) return -1;
172 return idx;
173 }
174
175 //__________________________________________________________________
176 Int_t BinCenter(UShort_t chip, UShort_t strip) const
177 {
178 return chip * 128 + strip;
179 }
180 //__________________________________________________________________
181 Int_t BinNumber(UShort_t chip, UShort_t strip) const
182 {
183 return 1 + BinCenter(chip, strip);
184 }
185 //__________________________________________________________________
186 Float_t GetSignal(UShort_t d, Char_t, UShort_t s, UShort_t t)
187 {
188 // Calculate signal value. Only the first 256 strips of FMD are
189 // filled.
190 if (d > 1 || s > 0 || t > 256) return 0;
191 Float_t c = 0.005;
192 Float_t x = (t % 128);
193 Float_t p = c * (64 * 64 - 128 * x + x * x) + 100;
194 Float_t n = (p-100)/40 + 1;
195
196 // In case we simulate a common mode noise component.
197 if (t < 128)
198 return fShift + gRandom->Gaus(p, n);
199
200 // In case all noise is uncorrelated.
201 return gRandom->Gaus(p, TMath::Sqrt(n*n + fCmn*fCmn));
202 }
203 //__________________________________________________________________
204 Bool_t ProcessRawDigit(AliFMDDigit* digit)
205 {
206 // From data
207 if (!digit) return kTRUE;
208
209 UShort_t d = digit->Detector();
210 Char_t r = digit->Ring();
211 UShort_t s = digit->Sector();
212 UShort_t t = digit->Strip();
213 return ProcessOne(d, r, s, t, digit->Counts());
214 }
215 //__________________________________________________________________
216 Bool_t ProcessUser(UShort_t d, Char_t r, UShort_t s, UShort_t t, Float_t v)
217 {
218 // From our simulator in GetValue
219 return ProcessOne(d, r, s, t, v);
220 }
221 //__________________________________________________________________
222 Bool_t ProcessOne(UShort_t d, Char_t r, UShort_t s, UShort_t t, Float_t v)
223 {
224 Int_t i = BinCenter(d, r, s, t);
225 if (i < 0) {
226 std::cout << "Invalid index " << i << " returned for FMD"
227 << d << r << '[' << s << ',' << t << ']' << std::endl;
228 return kFALSE;
229 }
230
231 fAdc->Fill(i, v);
232
233 return kTRUE;
234 }
235 //__________________________________________________________________
236 Bool_t End()
237 {
238 // At the end of each event, fill the summed histograms.
239 for (UShort_t chip = 0; chip < 400; chip++) {
240 Double_t chipAv = 0;
241 for (UShort_t strip = 0; strip < 128; strip++)
242 chipAv += fAdc->GetBinContent(1+BinCenter(chip, strip));
243 chipAv /= 128;
244
245 for (UShort_t strip = 0; strip < 128; strip++) {
246 Int_t i = BinCenter(chip, strip);
247 Double_t adc = fAdc->GetBinContent(i+1);
248 Double_t cAdc = chipAv - adc;
249 fSummedAdc->Fill(i, adc);
250 fSummedAdc2->Fill(i, adc * adc);
251 fSummedCenteredAdc->Fill(i, cAdc);
252 fSummedCenteredAdc2->Fill(i, cAdc*cAdc);
253 }
254 }
255
256 return AliFMDInput::End();
257 }
258 //__________________________________________________________________
259 Bool_t Finish()
260 {
261 if (fEventCount == 0) {
262 std::cerr << "No events!" << std::endl;
263 return kFALSE;
264 }
265
266 for (UInt_t bin = 1; bin <= 51200; bin++) {
267 if (bin % (51200 / 10) == 0)
268 std::cout << "Bin # " << bin << std::endl;
269 Double_t sumAdc = fSummedAdc->GetBinContent(bin);
270 Double_t sumAdc2 = fSummedAdc2->GetBinContent(bin);
271 Double_t avAdc = sumAdc / fEventCount;
272 Double_t avAdc2 = sumAdc2 / fEventCount;
273 Double_t noise2 = avAdc2 - avAdc * avAdc;
274
275 fPed->SetBinContent(bin, avAdc);
276 fNoise->SetBinContent(bin, noise2 >= 0 ? TMath::Sqrt(noise2) : 0);
277
278 Double_t sumCAdc = fSummedCenteredAdc->GetBinContent(bin);
279 Double_t sumCAdc2 = fSummedCenteredAdc2->GetBinContent(bin);
280 Double_t avCAdc = sumCAdc / fEventCount;
281 Double_t avCAdc2 = sumCAdc2 / fEventCount;
282 Double_t cNoise2 = avCAdc2 - avCAdc * avCAdc;
283 fCenteredPed->SetBinContent(bin, avCAdc);
284 fCorrectedNoise->SetBinContent(bin, TMath::Sqrt(cNoise2));
285 }
286
287 gStyle->SetPalette(1);
288 gStyle->SetOptTitle(0);
289 gStyle->SetCanvasColor(0);
290 gStyle->SetCanvasBorderSize(0);
291 gStyle->SetPadColor(0);
292 gStyle->SetPadBorderSize(0);
293
294 TCanvas* c = new TCanvas("c", "C", 800, 600);
295 c->SetFillColor(0);
296 c->SetBorderMode(0);
297 c->SetBorderSize(0);
298 c->SetTopMargin(0);
299 c->SetBottomMargin(0);
300 c->Divide(2, 1);
301
302 TPad* p1 = static_cast<TPad*>(c->cd(1));
303 p1->SetTopMargin(0.05);
304 p1->SetRightMargin(0.05);
305 if (!fReader) fNoise->GetXaxis()->SetRangeUser(-0.5, 255.5);
306 fNoise->SetMinimum(0);
307 fNoise->SetMaximum(1.5*fNoise->GetMaximum());
308 fNoise->DrawCopy();
309 fCorrectedNoise->SetMinimum(0);
310 fCorrectedNoise->SetFillColor(kBlue);
311 fCorrectedNoise->DrawCopy("same");
312 TLegend* l1 = new TLegend(0.2, 0.7, 0.945, 0.945);
313 l1->SetFillColor(0);
314 l1->SetFillStyle(0);
315 l1->SetBorderSize(0);
316 l1->AddEntry(fNoise, "Noise", "f");
317 l1->AddEntry(fCorrectedNoise, "Corrected noise", "f");
318 l1->Draw();
319
320
321 TPad* p2 = static_cast<TPad*>(c->cd(2));
322 p2->SetTopMargin(0.05);
323 p2->SetRightMargin(0.05);
324 if (!fReader) fPed->GetXaxis()->SetRangeUser(-0.5, 255.5);
325 fPed->SetMinimum(fCenteredPed->GetMinimum());
326 fPed->SetMaximum(1.5*fPed->GetMaximum());
327 fPed->DrawCopy();
328 fCenteredPed->SetFillColor(kBlue);
329 fCenteredPed->DrawCopy("same");
330 TLegend* l2 = new TLegend(0.2, 0.7, 0.945, 0.945);
331 l2->SetFillColor(0);
332 l2->SetFillStyle(0);
333 l2->SetBorderSize(0);
334 l2->AddEntry(fPed, "Pedestal", "f");
335 l2->AddEntry(fCenteredPed, "Pedestal - #bar{Pedestal}", "f");
336 l2->Draw();
337
338 c->Modified();
339 c->Update();
340 c->cd();
341 c->Print("cmn.png");
342
343 if (fOut) {
344 std::cout << "Flusing to disk ... " << std::flush;
345 fOut->cd();
346 fPed->Write(); //
347 fNoise->Write(); //
348 fCenteredPed->Write(); //
349 fCorrectedNoise->Write(); //
350 fSummedAdc->Write(); //
351 fSummedAdc2->Write(); //
352 fSummedCenteredAdc->Write(); //
353 fSummedCenteredAdc2->Write(); //
354 fAdc->Write(); //
355 fOut->Write();
356 fOut->Close();
357 std::cout << "done" << std::endl;
358 }
359 return kTRUE;
360 }
361
362 ClassDef(FindCommonModeNoise,0);
363};
364
365//____________________________________________________________________
366//
367// EOF
368//