Fixes for FMD DAs
[u/mrichter/AliRoot.git] / FMD / FMDutil / AliFMDPedestalDA.cxx
1 /**************************************************************************
2  * Copyright(c) 2004, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /** @file    AliMDPedestalDA.cxx
17     @author  Hans Hjersing Dalsgaard <canute@nbi.dk>
18     @date    Mon Mar 10 09:46:05 2008
19     @brief   Derived class for the pedestal detector algorithm.
20 */
21 //
22 // This class implements the virtual functions of the AliFMDBaseDA
23 // class.  The most important of these functions, FillChannels(..) and
24 // Analyse(...) collect and analyse the data of each channel. The
25 // resulting pedestal and noise values are written to a comma
26 // separated values (csv) file on the go. The csv files produced in
27 // this way are the basic input to the AliFMDPreprocessor.
28 //
29
30 #include "AliFMDPedestalDA.h"
31 #include "AliFMDAltroMapping.h"
32 #include "AliFMDParameters.h"
33 #include "AliFMDCalibPedestal.h"
34 #include "AliFMDDigit.h"
35 #include "AliLog.h"
36 #include <iostream>
37 #include <fstream>
38 #include <iomanip>
39 #include <TFile.h>
40 #include <TF1.h>
41 #include <TObject.h>
42 #include <TMath.h>
43 #include <TSystem.h>
44 #include <TDatime.h>
45 #include <TH2.h>
46 #include <TROOT.h>
47
48 //_____________________________________________________________________
49 ClassImp(AliFMDPedestalDA)
50
51 //_____________________________________________________________________
52 AliFMDPedestalDA::AliFMDPedestalDA() 
53   : AliFMDBaseDA(),
54     fCurrentChannel(1),
55     fPedSummary("PedestalSummary","pedestals",51200,0,51200),
56     fNoiseSummary("NoiseSummary","noise",51200,0,51200),
57     fZSfileFMD1(),
58     fZSfileFMD2(),
59     fZSfileFMD3(), 
60     fMinTimebin(3 * 4 * 3 * 16), // 3 ddls, 4 FECs, 3 Altros, 16 channels
61     fMaxTimebin(3 * 4 * 3 * 16), // 3 ddls, 4 FECs, 3 Altros, 16 channels
62     fSummaryFMD1i(0),
63     fSummaryFMD2i(0),
64     fSummaryFMD2o(0),
65     fSummaryFMD3i(0),
66     fSummaryFMD3o(0)
67 {
68   // Default constructor 
69   Rotate("peds.csv", 3);
70   fOutputFile.open("peds.csv");
71   Rotate("ddl3072.csv", 10);
72   fZSfileFMD1.open("ddl3072.csv");
73   Rotate("ddl3073.csv", 10);
74   fZSfileFMD2.open("ddl3073.csv");
75   Rotate("ddl3074.csv", 10);
76   fZSfileFMD3.open("ddl3074.csv");  
77   fDiagnosticsFilename = "diagnosticsPedestal.root";
78 }
79
80 //_____________________________________________________________________
81 AliFMDPedestalDA::AliFMDPedestalDA(const AliFMDPedestalDA & pedDA) : 
82   AliFMDBaseDA(pedDA),
83   fCurrentChannel(1),
84   fPedSummary("PedestalSummary","pedestals",51200,0,51200),
85   fNoiseSummary("NoiseSummary","noise",51200,0,51200),
86   fZSfileFMD1(),
87   fZSfileFMD2(),
88   fZSfileFMD3(),
89   fMinTimebin(pedDA.fMinTimebin),
90   fMaxTimebin(pedDA.fMaxTimebin),
91   fSummaryFMD1i(pedDA.fSummaryFMD1i),
92   fSummaryFMD2i(pedDA.fSummaryFMD2i),
93   fSummaryFMD2o(pedDA.fSummaryFMD2o),
94   fSummaryFMD3i(pedDA.fSummaryFMD3i),
95   fSummaryFMD3o(pedDA.fSummaryFMD3o)
96 {
97   // Copy constructor 
98 }
99
100 //_____________________________________________________________________
101 AliFMDPedestalDA::~AliFMDPedestalDA() 
102
103   // Destructor.
104 }
105
106 //_____________________________________________________________________
107 void AliFMDPedestalDA::Init() 
108
109   // Initialise 
110   SetRequiredEvents(1000);
111   fMinTimebin.Reset(1024);
112   fMaxTimebin.Reset(-1);
113
114
115 }
116
117 //_____________________________________________________________________
118 void AliFMDPedestalDA::AddChannelContainer(Array* sampleArray, 
119                                            UShort_t det, 
120                                            Char_t   ring, 
121                                            UShort_t sec, 
122                                            UShort_t strip) 
123 {
124   // Add a channel to the containers.
125   //
126   // Parameters:
127   //     sectorArray  Array of sectors
128   //     det          Detector 
129   //     ring         Ring
130   //     sec          Sector 
131   //     strip        Strip
132   AliFMDParameters* pars        = AliFMDParameters::Instance();
133   UInt_t            samples     = pars->GetSampleRate(det, ring, sec, strip);
134   for (UInt_t sample = 0; sample < samples; sample++) {
135     TString name(Form("FMD%d%c[%02d,%03d]_%d", det,ring,sec,strip,sample));
136     TH1S* hSample = new TH1S(name.Data(),name.Data(), 1024,-.5,1023.5);
137     hSample->SetXTitle("ADC");
138     hSample->SetYTitle("Events");
139     hSample->SetDirectory(0);
140     hSample->ResetBit(TObject::kMustCleanup);
141     sampleArray->AddAtAndExpand(hSample, sample);
142   }
143 }
144
145 //_____________________________________________________________________
146 void AliFMDPedestalDA::AddSectorSummary(Array* sectorArray, 
147                                         UShort_t det, 
148                                         Char_t   ring, 
149                                         UShort_t sec, 
150                                         UShort_t nStr) 
151 {
152   TH1F* sumPed = new TH1F("Pedestals", 
153                           Form("Summary of pedestals in FMD%d%c[%02d]", 
154                                det, ring, sec), 
155                           nStr, -.5, nStr-.5);
156   sumPed->SetXTitle("Strip");
157   sumPed->SetYTitle("Pedestal [ADC]");
158   sumPed->SetDirectory(0);
159   sumPed->ResetBit(TObject::kMustCleanup);
160   
161   TH1F* sumNoise = static_cast<TH1F*>(sumPed->Clone("Noise"));
162   sumNoise->SetYTitle("Noise [ADC]");
163   sumNoise->SetDirectory(0);
164   sumNoise->ResetBit(TObject::kMustCleanup);
165   
166   Int_t n = sectorArray->GetEntriesFast();
167   sectorArray->AddAtAndExpand(sumPed,   n + kPedestalOffset - 1);
168   sectorArray->AddAtAndExpand(sumNoise, n + kNoiseOffset - 1);
169 }
170
171 //_____________________________________________________________________
172 void AliFMDPedestalDA::FillChannels(AliFMDDigit* digit) 
173 {
174   // Fill ADC values from a digit into the corresponding histogram.
175   //
176   // Parameters:
177   //    digit Digit to fill ADC values for.
178   UShort_t det   = digit->Detector();
179   Char_t   ring  = digit->Ring();
180   UShort_t sec   = digit->Sector();
181   UShort_t strip = digit->Strip();
182   
183   AliFMDParameters* pars     = AliFMDParameters::Instance();
184   UInt_t            samples  = pars->GetSampleRate(det, ring, sec, strip);
185   for (UInt_t sample = 0; sample < samples; sample++) {
186     TH1S* hSample = GetChannel(det, ring, sec, strip, sample);
187     if (!hSample) continue;
188     
189     hSample->Fill(digit->Count(sample));
190   }
191   
192 }
193
194 //_____________________________________________________________________
195 void AliFMDPedestalDA::MakeSummary(UShort_t det, Char_t ring)
196 {
197   //Create summary hists for FMD pedestals
198   // std::cout << "Making summary for FMD" << det << ring << " ..." 
199   //           << std::endl;
200   switch (det) { 
201   case 1: 
202     fSummaryFMD1i = MakeSummaryHistogram("ped", "Pedestals", det, ring);
203     break;
204   case 2:
205     switch (ring) { 
206     case 'I': case 'i':
207       fSummaryFMD2i = MakeSummaryHistogram("ped", "Pedestals", det, ring);
208       break;
209     case 'O': case 'o':
210       fSummaryFMD2o = MakeSummaryHistogram("ped", "Pedestals", det, ring);
211       break;
212     }
213     break;
214   case 3:
215     switch (ring) { 
216     case 'I': case 'i':
217       fSummaryFMD3i = MakeSummaryHistogram("ped", "Pedestals", det, ring);
218       break;
219     case 'O': case 'o':
220       fSummaryFMD3o = MakeSummaryHistogram("ped", "Pedestals", det, ring);
221       break;
222     }
223     break;
224   }
225 }
226
227 //_____________________________________________________________________
228 void AliFMDPedestalDA::Analyse(UShort_t det, 
229                                Char_t   ring, 
230                                UShort_t sec, 
231                                UShort_t strip) 
232 {
233   // Analyse a strip.  That is, compute the mean and spread of the ADC
234   // spectra for all strips.  Also output on files the values. 
235   //
236   // Parameters:
237   //     det   Detector
238   //     ring  Ring
239   //     sec   Sector 
240   //     strip Strip.
241   AliFMDParameters* pars     = AliFMDParameters::Instance();
242   TH2* summary = 0;
243   switch (det) { 
244   case 1: summary = fSummaryFMD1i; break;
245   case 2: 
246     switch (ring) { 
247     case 'I':  summary = fSummaryFMD2i; break;
248     case 'O':  summary = fSummaryFMD2o; break;
249     }
250     break;
251   case 3:
252     switch (ring) { 
253     case 'I':  summary = fSummaryFMD3i; break;
254     case 'O':  summary = fSummaryFMD3o; break;
255     }
256     break;
257   }
258
259   UInt_t samples  = pars->GetSampleRate(det, ring, sec, strip);
260   for (UShort_t sample = 0; sample < samples; sample++) {
261     TH1S* hChannel = GetChannel(det, ring, sec, strip,sample);
262     if(!hChannel || hChannel->GetEntries() == 0) {
263       //AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
264       //            det,ring,sec,strip));
265       return;
266     }
267     
268     AliDebugF(50, "Fitting FMD%d%c[%02d,%03d] with %d entries",
269               det,ring,sec,strip, int(hChannel->GetEntries()));
270     TF1* fitFunc = new TF1("fitFunc","gausn",0,300);
271     fitFunc->ResetBit(TObject::kMustCleanup);
272     fitFunc->SetParameters(100,100,1);
273     hChannel->Fit(fitFunc,"Q0","",10,200);
274     hChannel->GetListOfFunctions()->Remove(fitFunc);
275     gROOT->GetListOfFunctions()->Remove(fitFunc);
276
277     Float_t mean = hChannel->GetMean();
278     Float_t rms  = hChannel->GetRMS();
279     
280     hChannel->GetXaxis()->SetRangeUser(mean-5*rms,mean+5*rms);
281     mean = hChannel->GetMean();
282     rms  = hChannel->GetRMS();
283   
284     
285     UShort_t ddl, board, altro, channel;
286     UShort_t timebin;
287     
288     pars->Detector2Hardware(det,ring,sec,strip,sample,
289                             ddl,board,altro,channel,timebin);
290     Int_t idx = HWIndex(ddl, board, altro, channel);
291     if (idx >= 0) { 
292       fMinTimebin[idx] = TMath::Min(Short_t(timebin),   fMinTimebin[idx]);
293       fMaxTimebin[idx] = TMath::Max(Short_t(timebin+1), fMaxTimebin[idx]);
294     }
295     
296     std::ostream* zsFile = 0;
297     switch(det) {
298     case 1:  zsFile = &fZSfileFMD1; break;
299     case 2:  zsFile = &fZSfileFMD2; break;
300     case 3:  zsFile = &fZSfileFMD3; break;
301     default: AliWarning("Unknown sample!"); break;
302       
303     }
304     *zsFile << board   << ',' 
305             << altro   << ',' 
306             << channel << ',' 
307             << timebin << ','
308             << mean    << ',' 
309             << rms     << "\n";
310     
311     Float_t chi2ndf = 0;
312     
313     
314     if(fitFunc->GetNDF())
315       chi2ndf = fitFunc->GetChisquare() / fitFunc->GetNDF();
316     
317     
318     Int_t sampleToWrite = 2;
319     if      (samples == 2) sampleToWrite = 1;
320     else if (samples <  2) sampleToWrite = 0;
321     
322     hChannel->GetXaxis()->SetRange(1,1024);
323
324     if(sample != sampleToWrite) continue;
325     
326     
327     fOutputFile << det                         << ','
328                 << ring                        << ','
329                 << sec                         << ','
330                 << strip                       << ','
331                 << mean                        << ','
332                 << rms                         << ','
333                 << fitFunc->GetParameter(1)    << ','
334                 << fitFunc->GetParameter(2)    << ','
335                 << chi2ndf                     <<"\n";
336
337     delete fitFunc;
338
339     if (summary) { 
340       Int_t bin = summary->FindBin(sec, strip);
341       summary->SetBinContent(bin, mean);
342       summary->SetBinError(bin, rms);
343     }
344
345     if(fSaveHistograms  ) {
346       TH1F* sumPed   = GetSectorSummary(det, ring, sec, true);
347       TH1F* sumNoise = GetSectorSummary(det, ring, sec, false);
348       sumPed->SetBinContent(strip+1, mean);
349       sumPed->SetBinError(strip+1, rms);
350       sumNoise->SetBinContent(strip+1, rms);
351
352       fPedSummary.SetBinContent(fCurrentChannel,mean);
353       fNoiseSummary.SetBinContent(fCurrentChannel,rms);
354       fCurrentChannel++;
355     }
356   }
357 }
358
359 //_____________________________________________________________________
360 void AliFMDPedestalDA::Terminate(TFile* diagFile) 
361 {
362   // Called at the end of a job.  Fills in missing time-bins and
363   // closes output files
364   if(fSaveHistograms && diagFile) {
365     diagFile->cd();
366     
367     fPedSummary.Write();
368     fNoiseSummary.Write();
369   }
370   AliFMDAltroMapping* map = AliFMDParameters::Instance()->GetAltroMap();
371   for (Int_t i = 0; i < 3; i++) { 
372     std::ofstream& out = (i == 0 ? fZSfileFMD1 : 
373                           i == 1 ? fZSfileFMD2 : 
374                           fZSfileFMD3);
375     if (out.is_open() && fSeenDetectors[i]) { 
376       FillinTimebins(out, map->Detector2DDL(i+1));
377     }
378     if (!fSeenDetectors[i]) {
379       TString n(Form("ddl%d.csv",3072+map->Detector2DDL(i+1)));
380       gSystem->Unlink(n.Data());
381     }
382   }
383   
384 }
385
386 //_____________________________________________________________________
387 void AliFMDPedestalDA::FillinTimebins(std::ofstream& out, UShort_t /*ddl*/)
388 {
389   // 
390   // Fill missing timebins
391   //
392 #if 0
393   unsigned short  boards[] = { 0x0, 0x1, 0x10, 0x11, 0xFFFF };
394   unsigned short* board    = boards;
395   while ((*boards) != 0xFFFF) { 
396     for (UShort_t altro = 0; altro < 3; altro++) { 
397       for (UShort_t channel = 0; channel < 16; channel++) { 
398         Int_t idx = HWIndex(ddl, *board, altro, channel);
399         if (idx < 0) { 
400           AliWarningF("Invalid index for %4d/0x%02x/0x%x/0x%x: %d", 
401                       ddl, *board, altro, channel, idx);
402           continue;
403         }
404         Short_t min = fMinTimebin[idx];
405         Short_t max = fMaxTimebin[idx];
406         
407         // Channel not seen at all. 
408         if (min > 1023 || max < 0) continue;
409
410         out << "# Extra timebins for 0x" << std::hex 
411             << board << ',' << altro << ',' << channel 
412             << " got time-bins " << min << " to " << max-1 
413             << std::dec << std::endl;
414         
415         for (UShort_t t = 15; t < min; t++) 
416           // Write a phony line 
417           out << board << "," << altro << "," << channel << "," 
418               << t << "," << 1023 << "," << 0 << std::endl;
419
420         for (UShort_t t = max; t < 1024; t++) 
421           // Write a phony line 
422           out << board << "," << altro << "," << channel << "," 
423               << t << "," << 1023 << "," << 0 << std::endl;
424       } // channel loop
425     } // altro loop
426   } // board loop
427   // Write trailer, and close 
428 #endif
429   out.write("# EOF\n", 6);
430   out.close();
431 }
432
433 //_____________________________________________________________________
434 void AliFMDPedestalDA::WriteHeaderToFile() 
435 {
436   //
437   // Write headers to output files 
438   //
439   AliFMDParameters* pars       = AliFMDParameters::Instance();
440   fOutputFile.write(Form("# %s \n",pars->GetPedestalShuttleID()),13);
441   TDatime now;
442   fOutputFile << "# This file created from run # " << fRunno 
443               << " @ " << now.AsString() << std::endl;
444   fOutputFile.write("# Detector, "
445                     "Ring, "
446                     "Sector, "
447                     "Strip, "
448                     "Pedestal, "
449                     "Noise, "
450                     "Mu, "
451                     "Sigma, "
452                     "Chi2/NDF \n", 71);
453
454   std::ostream* zss[] = { &fZSfileFMD1, &fZSfileFMD2, &fZSfileFMD3, 0 };
455   for (size_t i = 0; i < 3; i++)  {
456     *(zss[i]) << "# FMD " << (i+1) << " pedestals \n"
457               << "# board, "
458               << "altro, "
459               << "channel, "
460               << "timebin, "
461               << "pedestal, "
462               << "noise\n";
463     *(zss[i]) << "# This file created from run # " << fRunno 
464               << " @ " << now.AsString() << std::endl;
465   }
466 }
467
468 //_____________________________________________________________________
469 TH1S* AliFMDPedestalDA::GetChannel(UShort_t det, 
470                                    Char_t   ring, 
471                                    UShort_t sec, 
472                                    UShort_t strip,
473                                    UInt_t   sample) 
474 {
475   // Get the histogram corresponding to a strip sample.
476   // 
477   // Parameters:
478   //     det    Detector
479   //     ring   Ring
480   //     sec    Sector
481   //     strip  Strip
482   //     sample Sample
483   //
484   // Return:
485   //     ADC spectra of a strip.
486   Array* sampleArray = GetStripArray(det, ring, sec, strip);
487   if (!sampleArray) return 0;
488   TH1S*      hSample     = static_cast<TH1S*>(sampleArray->At(sample));
489   if (!hSample) {
490     AliErrorF("No channel histogram for FMD%d%c[%02d,%03d]_%d", 
491               det, ring, sec, strip, sample);
492     sampleArray->ls();
493     AliErrorF("Path is %s <- %s <- %s <- %s", 
494               sampleArray->GetName(),
495               GetSectorArray(det, ring, sec)->GetName(), 
496               GetRingArray(det, ring)->GetName(),
497               GetDetectorArray(det)->GetName());
498     
499   }
500   return hSample;
501   
502 }
503 //_____________________________________________________________________
504 TH1F* AliFMDPedestalDA::GetSectorSummary(UShort_t det, 
505                                          Char_t   ring, 
506                                          UShort_t sec, 
507                                          Bool_t   pedNotNoise) 
508 {
509   Array*     secArray    = GetSectorArray(det, ring, sec);
510   Int_t      n           = secArray->GetEntriesFast();
511   Int_t      i           = n - (pedNotNoise ? kNoiseOffset : kPedestalOffset);
512   return static_cast<TH1F*>(secArray->At(i));
513 }
514
515 //_____________________________________________________________________
516 //
517 //EOF
518 //