bug fix
[u/mrichter/AliRoot.git] / FMD / scripts / DrawCalibRaw.C
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 <TH1D.h>
15 #include <AliCDBManager.h>
16 #include <AliFMDHit.h>
17 #include <AliFMDDigit.h>
18 #include <AliFMDInput.h>
19 #include <AliFMDEdepMap.h>
20 #include <AliFMDParameters.h>
21 #include <AliFMDCalibStripRange.h>
22 #include <AliFMDCalibSampleRate.h>
23 #include <AliFMDCalibGain.h>
24 #include <AliFMDCalibPedestal.h>
25 #include <AliFMDRawReader.h>
26 #include <iostream>
27 #include <fstream>
28 #include <TStyle.h>
29 #include <TArrayF.h>
30 #include <AliLog.h>
31 #include <TSystem.h>
32 #include <TFile.h>
33
34 /** @class DrawDigits
35     @brief Draw hit energy loss versus digit ADC
36     @code 
37     Root> .L Compile.C
38     Root> Compile("DrawCalibRaw.C")
39     Root> DrawCalibRaw cr;
40     Root> cr.SetRawFile("file");
41     Root> cr.Run();
42     @endcode
43     @ingroup FMD_script
44  */
45 class DrawCalibRaw : public AliFMDInput
46 {
47 private:
48   Double_t   fFactor;
49   TString    fCalibDir;
50   TH1D*      fELoss; // Histogram 
51   TObjArray* fDets;
52   TFile*     fOut;
53   Bool_t     fHasData;
54   Int_t      fGotNEvents;
55
56 public:
57   //__________________________________________________________________
58   DrawCalibRaw(const       char*  file, 
59                const char* calibDir    = 0,
60                Double_t    noiseFactor = 5, 
61                Bool_t      save        = kTRUE,
62                Int_t       m           = 420, 
63                Double_t    mmin        = -0.5, 
64                Double_t    mmax        = 20.5) 
65     : AliFMDInput(), 
66       fFactor(noiseFactor),
67       fCalibDir(""), 
68       fELoss(0), 
69       fDets(0), 
70       fOut(0), 
71       fHasData(kFALSE),
72       fGotNEvents(0)
73   { 
74     if (calibDir) fCalibDir = calibDir;
75     AddLoad(kRawCalib);
76     fELoss = new TH1D("eLoss", "Scaled Energy loss", m, mmin, mmax);
77     fELoss->SetXTitle("#Delta E/#Delta E_{mip}");
78     fELoss->Sumw2();
79     fELoss->SetDirectory(0);
80     
81     SetRawFile(file);
82
83     if (!save) return;
84     fDets = new TObjArray(3);
85
86   }
87   //__________________________________________________________________
88   Bool_t CheckFile(const char* prefix, int number, TString& f)
89   {
90     f = (Form("%s%d.csv", prefix, number));
91     f = gSystem->Which(fCalibDir.Data(), f.Data());
92     return !f.IsNull();
93   }
94
95   //__________________________________________________________________
96   Bool_t Init()
97   {
98     AliCDBManager* cdb = AliCDBManager::Instance();
99     cdb->SetRun(0);
100     cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
101
102     AliFMDCalibStripRange* range = 0;
103     AliFMDCalibSampleRate* rate  = 0;
104     AliFMDCalibPedestal*   peds  = 0;
105     AliFMDCalibGain*       gains = 0;
106     for (Int_t i = 1; i <= 3; i++) { 
107       TString f;
108       if (CheckFile("conditions", i, f)) {
109         Info("Init", "Reading conditions for FMD%d from %s", i, f.Data());
110         if (!range) range = new AliFMDCalibStripRange;
111         if (!rate)  rate  = new AliFMDCalibSampleRate;
112         std::ifstream in(f.Data());
113         range->ReadFromFile(in);
114         rate->ReadFromFile(in);
115       }
116       if (CheckFile("peds", i, f)) {
117         Info("Init", "Reading pedestals for FMD%d from %s", i, f.Data());
118         if (!peds) peds = new AliFMDCalibPedestal;
119         std::ifstream in(f.Data());
120         peds->ReadFromFile(in);
121       }
122       if (CheckFile("gains", i, f)) {
123         Info("Init", "Reading gains for FMD%d from %s", i, f.Data());
124         if (!gains) gains = new AliFMDCalibGain;
125         std::ifstream in(f.Data());
126         gains->ReadFromFile(in);
127       }
128     }
129
130     Int_t mask = (AliFMDParameters::kDeadMap|
131                   AliFMDParameters::kZeroSuppression|
132                   AliFMDParameters::kAltroMap);
133
134     if (!range) mask |= AliFMDParameters::kStripRange;
135     if (!rate)  mask |= AliFMDParameters::kSampleRate;
136     if (!peds)  mask |= AliFMDParameters::kPedestal;
137     if (!gains) mask |= AliFMDParameters::kPulseGain;
138
139     AliFMDParameters* pars = AliFMDParameters::Instance();
140     pars->Init(kFALSE, mask);
141
142     if (range)  pars->SetStripRange(range);
143     if (rate)   pars->SetSampleRate(rate);
144     if (peds)   pars->SetPedestal(peds);
145     if (gains)  pars->SetGain(gains);
146     
147     Bool_t ret = AliFMDInput::Init();
148     
149     if (!fDets) return ret;
150
151     fOut  = TFile::Open(Form("histo_%s", fRawFile.Data()), "RECREATE");
152
153     Int_t m    = fELoss->GetXaxis()->GetNbins();
154     Int_t mmin = fELoss->GetXaxis()->GetXmin();
155     Int_t mmax = fELoss->GetXaxis()->GetXmax();
156     for (Int_t d = 1; d <= 3; d++) {
157       Int_t      nRng = (d == 1 ? 1 : 2);
158       TObjArray* det  = 0;
159       fDets->AddAt(det = new TObjArray(nRng), d-1);
160       det->SetName(Form("FMD%d", d));
161       TDirectory* detD = fOut->mkdir(det->GetName());
162       for (Int_t q = 0; q < nRng; q++) { 
163         Char_t r       = q == 0 ? 'I' : 'O';
164         Int_t  nSec    = q == 0 ?  20 :  40;
165         Int_t  nStr    = q == 0 ? 512 : 256;
166         TObjArray* rng = 0;
167         det->AddAt(rng = new TObjArray(nSec), q);
168         rng->SetName(Form("FMD%d%c", d, r));
169         TDirectory* rngD = detD->mkdir(rng->GetName());
170         for (Int_t s = 0; s < nSec; s++) { 
171           TObjArray* sec = 0;
172           rng->AddAt(sec = new TObjArray(nStr), s);
173           sec->SetName(Form("FMD%d%c_%02d", d, r, s));
174           TDirectory* secD = rngD->mkdir(sec->GetName());
175           for (Int_t t = 0; t < nStr; t++) { 
176             secD->cd();
177             TH1* str = new TH1D(Form("FMD%d%c_%02d_%03d", d, r, s, t), 
178                                  Form("Scaled energy loss in FMD%d%c[%2d,%3d]",
179                                       d, r, s, t), m, mmin, mmax);
180             str->SetXTitle("#Delta E/#Delta E_{mip}");
181             // str->SetDirectory(secD);
182             sec->AddAt(str, t);
183           }
184         }
185       }
186     }
187
188     return ret;
189   }
190
191   //__________________________________________________________________
192   Bool_t Begin(Int_t e)
193   {
194     fHasData = kFALSE;
195     return AliFMDInput::Begin(e);
196   }
197
198   //__________________________________________________________________
199   Bool_t ProcessRawCalibDigit(AliFMDDigit* digit)
200   {
201     if (!digit) return kTRUE;
202
203     AliFMDParameters* parm = AliFMDParameters::Instance();
204     UShort_t d             =  digit->Detector();
205     Char_t   r             =  digit->Ring();
206     UShort_t s             =  digit->Sector();
207     UShort_t t             =  digit->Strip();
208     Double_t gain          =  parm->GetPulseGain(d, r, s, t);
209     Double_t ped           =  parm->GetPedestal(d, r, s, t);
210     Double_t pedW          =  parm->GetPedestalWidth(d, r, s, t);
211     Double_t adc           =  digit->Counts();
212     Double_t threshold     =  pedW * fFactor;
213     if (gain < 0.1 || gain > 10) return kTRUE;
214     if (pedW > 10) { 
215       Warning("ProcessRawCalibDigit", "FMD%d%c[%2d,%3d] is noisy: %f",
216               d, r, s, t, pedW);
217       return kTRUE;
218     }
219
220     if (fFMDReader && fFMDReader->IsZeroSuppressed(d-1))
221       adc += fFMDReader->NoiseFactor(d-1) * pedW;
222     else 
223       threshold            += ped;
224
225     if (adc < threshold) return kTRUE;
226     
227     Double_t mult = (adc-ped) / (gain * parm->GetDACPerMIP());
228
229     fHasData = kTRUE;
230     fELoss->Fill(mult);
231
232     // if (t >= 10) return kTRUE;
233     TObjArray* det = static_cast<TObjArray*>(fDets->At(d-1));
234     TObjArray* rng = static_cast<TObjArray*>(det->At(r == 'I' ? 0 : 1));
235     TObjArray* sec = static_cast<TObjArray*>(rng->At(s));
236     TH1*       str = static_cast<TH1*>(sec->At(t));
237     str->Fill(mult);
238
239     return kTRUE;
240   }
241   //__________________________________________________________________
242   Bool_t End()
243   {
244     if (fHasData) fGotNEvents++;
245     return AliFMDInput::End();
246   }
247   //__________________________________________________________________
248   Bool_t Finish()
249   {
250     std::cout << "A total of " << fGotNEvents << " with data" << std::endl;
251     gStyle->SetPalette(1);
252     gStyle->SetOptTitle(0);
253     gStyle->SetCanvasColor(0);
254     gStyle->SetCanvasBorderSize(0);
255     gStyle->SetPadColor(0);
256     gStyle->SetPadBorderSize(0);
257     fELoss->SetStats(kFALSE);
258     fELoss->SetFillColor(kRed);
259     fELoss->SetFillStyle(3001);
260     fELoss->Scale(1. / fELoss->GetEntries());
261     fELoss->DrawCopy("e1 bar");
262
263     if (fDets && fOut) { 
264       std::cout << "Flusing to disk ... " << std::flush;
265       fOut->cd();
266       fELoss->Write();
267       fOut->Write();
268       fOut->Close();
269       std::cout << "done" << std::endl;
270     }
271     return kTRUE;
272   }
273
274   ClassDef(DrawCalibRaw,0);
275 };
276
277 //____________________________________________________________________
278 //
279 // EOF
280 //