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