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