Fixes for FMD DAs
[u/mrichter/AliRoot.git] / FMD / FMDutil / AliFMDGainDA.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    AliFMDGainDA.cxx
17     @author  Hans Hjersing Dalsgaard <canute@nbi.dk>
18     @date    Mon Mar 13 13:46:05 2008
19     @brief   Derived class for the pulse gain detector algorithm.
20 */
21 //____________________________________________________________________
22 //
23 // This class contains the implementation of the gain detector
24 // algorithms (DA) for the FMD.  The data is collected in histograms
25 // that are reset for each pulse length after the mean and standard
26 // deviation are put into a TGraphErrors object. After a certain
27 // number of pulses (usually 8) the graph is fitted to a straight
28 // line. The gain is then slope of this line as it combines the known
29 // pulse and the response of the detector.
30 //
31 #include "AliFMDGainDA.h"
32 #include "AliFMDAltroMapping.h"
33 #include "AliFMDParameters.h"
34 #include "AliFMDCalibGain.h"
35 #include "AliFMDDigit.h"
36 #include "AliLog.h"
37 #include <TFile.h>
38 #include <TF1.h>
39 #include <TH1S.h>
40 #include <TGraphErrors.h>
41 #include <TDatime.h>
42 #include <TH2.h>
43 #include <TROOT.h>
44
45 //_____________________________________________________________________
46 ClassImp(AliFMDGainDA)
47 #if 0 // Do not delete - here to let Emacs indent properly
48 ;
49 #endif
50
51 //_____________________________________________________________________
52 AliFMDGainDA::AliFMDGainDA() 
53   : AliFMDBaseDA(),
54     fHighPulse(256), 
55     fEventsPerChannel(10),
56     fCurrentPulse(10),
57     fCurrentChannel(10),
58     fNumberOfStripsPerChip(128),
59     fSummaryGains("GainsSummary","Summary of gains",51200,0,51200),
60     fCurrentSummaryStrip(1),
61     fGainFMD1i(0),
62     fGainFMD2i(0),
63     fGainFMD2o(0),
64     fGainFMD3i(0),
65     fGainFMD3o(0),
66     fChi2FMD1i(0),
67     fChi2FMD2i(0),
68     fChi2FMD2o(0),
69     fChi2FMD3i(0),
70     fChi2FMD3o(0)
71 {
72   // Constructor 
73   // 
74   // Parameters: 
75   //   None
76   fCurrentPulse.Reset(0);
77   fCurrentChannel.Reset(0);
78   fOutputFile.open("gains.csv");
79   fDiagnosticsFilename = "diagnosticsGain.root";
80 }
81
82 //_____________________________________________________________________
83 AliFMDGainDA::AliFMDGainDA(const AliFMDGainDA & gainDA) 
84   : AliFMDBaseDA(gainDA),
85     fHighPulse(gainDA.fHighPulse),
86     fEventsPerChannel(gainDA.fEventsPerChannel),
87     fCurrentPulse(gainDA.fCurrentPulse),
88     fCurrentChannel(gainDA.fCurrentChannel),
89     fNumberOfStripsPerChip(gainDA.fNumberOfStripsPerChip),
90     fSummaryGains(gainDA.fSummaryGains),
91     fCurrentSummaryStrip(gainDA.fCurrentSummaryStrip),
92     fGainFMD1i(gainDA.fGainFMD1i),
93     fGainFMD2i(gainDA.fGainFMD2i),
94     fGainFMD2o(gainDA.fGainFMD2o),
95     fGainFMD3i(gainDA.fGainFMD3i),
96     fGainFMD3o(gainDA.fGainFMD3o),
97     fChi2FMD1i(gainDA.fChi2FMD1i),
98     fChi2FMD2i(gainDA.fChi2FMD2i),
99     fChi2FMD2o(gainDA.fChi2FMD2o),
100     fChi2FMD3i(gainDA.fChi2FMD3i),
101     fChi2FMD3o(gainDA.fChi2FMD3o)
102 {  
103   // Copy Constructor 
104   // 
105   // Parameters: 
106   //   gainDA   Object to copy from 
107   fCurrentPulse.Reset(0);
108   fCurrentChannel.Reset(0);
109 }
110
111 //_____________________________________________________________________
112 AliFMDGainDA::~AliFMDGainDA() 
113 {
114   // Destructor 
115   // 
116   // Parameters: 
117   //   None
118 }
119
120 //_____________________________________________________________________
121 void AliFMDGainDA::Init() 
122 {
123   // Initialize 
124   // 
125   // Parameters: 
126   //   None
127   Int_t nEventsRequired = 0;
128   
129   for(Int_t idx = 0;idx<fEventsPerChannel.GetSize();idx++) {
130     Int_t nEvents = 0;
131     if(fPulseSize.At(idx))
132       nEvents = (fPulseLength.At(idx)*fHighPulse) / fPulseSize.At(idx);
133     fEventsPerChannel.AddAt(nEvents,idx);
134     if(nEvents>nEventsRequired) 
135       nEventsRequired = nEvents * fNumberOfStripsPerChip;    
136   }
137   SetRequiredEvents(nEventsRequired); 
138 }
139
140 //_____________________________________________________________________
141 void AliFMDGainDA::InitContainer(TDirectory* dir) 
142 {
143   AliFMDBaseDA::InitContainer(dir);
144   
145   for(UShort_t det=1;det<=3;det++) {
146     UShort_t sr = (det == 1 ? 1 : 0);
147     for (UShort_t ir = sr; ir < 2; ir++) {
148       Char_t   ring = (ir == 0 ? 'O' : 'I');
149       UShort_t nsec = (ir == 0 ? 40  : 20);
150       UShort_t nva  = (ir == 0 ? 2 : 4);
151       for(UShort_t sec =0; sec < nsec;  sec++)  {
152         Array* sectorArray = GetSectorArray(det, ring, sec);
153         Array* cache = new Array(nva);
154         cache->SetName("Cache");
155         cache->SetOwner();
156         Int_t n = sectorArray->GetEntriesFast();
157         sectorArray->AddAtAndExpand(cache, n);
158         for(UShort_t va = 0; va < nva; va++) {
159           TH1S* hChannel = new TH1S(Form("FMD%d%c[%02d]_va%d",
160                                          det,ring,sec,va),
161                                     Form("FMD%d%c[%02d] VA%d Cache",
162                                          det,ring,sec,va),
163                                     1024,-.5,1023.5);
164           hChannel->SetDirectory(0);
165           cache->AddAtAndExpand(hChannel,va);
166         }
167       }
168     }
169   }
170 }
171
172 //_____________________________________________________________________
173 void AliFMDGainDA::AddChannelContainer(Array* stripArray, 
174                                        UShort_t det  , 
175                                        Char_t   ring,  
176                                        UShort_t sec, 
177                                        UShort_t strip) 
178 {  
179   // Make a channel container 
180   // 
181   // Parameters: 
182   //  sectorArray    Sectors 
183   //  det            Detector number
184   //  ring           Ring identifier 
185   //  sec            Sector number
186   //  strip          Strip number
187
188   AliFMDParameters* pars     = AliFMDParameters::Instance();
189   UShort_t          board    = pars->GetAltroMap()->Sector2Board(ring, sec);
190   Int_t             halfring = GetHalfringIndex(det,ring,board/16);
191   Int_t             dPulse   = fPulseSize.At(halfring);
192   Int_t             nPulses  = dPulse > 0 ? fHighPulse / dPulse : 0;
193
194   TGraphErrors* gChannel  = new TGraphErrors(nPulses);
195   gChannel->SetName(Form("FMD%d%c[%02d,%03d]", det, ring, sec, strip));
196   gChannel->SetTitle(Form("FMD%d%c[%02d,%03d] ADC vs DAC", 
197                           det, ring, sec, strip));
198   stripArray->AddAtAndExpand(gChannel,0);
199 }
200
201 //_____________________________________________________________________
202 void AliFMDGainDA::AddSectorSummary(Array* sectorArray, 
203                                     UShort_t det, 
204                                     Char_t   ring, 
205                                     UShort_t sec, 
206                                     UShort_t nStr) 
207 {
208   TH1F* summary = new TH1F("Summary", Form("Summary of gains in FMD%d%c[%02d]", 
209                                            det, ring, sec), 
210                            nStr, -.5, nStr-.5);
211   summary->SetXTitle("Strip");
212   summary->SetYTitle("Gain [ADC/DAC]");
213   summary->SetDirectory(0);
214
215   Int_t n = sectorArray->GetEntriesFast();
216   sectorArray->AddAtAndExpand(summary, n);
217 }
218
219 //_____________________________________________________________________
220 void AliFMDGainDA::FillChannels(AliFMDDigit* digit) 
221 {
222   // Fill data into histogram
223   // 
224   // Parameters: 
225   //   digit  Digit to get the data from
226
227   UShort_t det   = digit->Detector();
228   Char_t   ring  = digit->Ring();
229   UShort_t sec   = digit->Sector();
230   UShort_t strip = digit->Strip();
231   
232   //Strip is always seen as the first in a VA chip. All other strips are junk.
233   //Strips are counted from zero on even sectors and from 511 on odd sectors...
234    
235   if((sec%2)     && ((strip+1) % fNumberOfStripsPerChip)) return;
236   if(((sec+1)%2) && (strip % fNumberOfStripsPerChip)) return;
237   
238   UShort_t vaChip   = strip / fNumberOfStripsPerChip; 
239   TH1S*    hChannel = GetChannelHistogram(det, ring, sec, vaChip);
240   hChannel->Fill(digit->Counts());
241   UpdatePulseAndADC(det,ring,sec,strip);
242 }
243
244 //_____________________________________________________________________
245 void AliFMDGainDA::MakeSummary(UShort_t det, Char_t ring)
246 {
247   // 
248   // Create summary hists for FMD gains and chi2 of the fits
249   //
250   switch (det) { 
251   case 1: 
252     fGainFMD1i = MakeSummaryHistogram("gain", "Gains", det, ring);
253     fChi2FMD1i = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
254     break;
255   case 2:
256     switch (ring) { 
257     case 'I': case 'i':
258       fGainFMD2i = MakeSummaryHistogram("gain", "Gains", det, ring);
259       fChi2FMD2i = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
260       break;
261     case 'O': case 'o':
262       fGainFMD2o = MakeSummaryHistogram("gain", "Gains", det, ring);
263       fChi2FMD2o = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
264       break;
265     }
266     break;
267   case 3:
268     switch (ring) { 
269     case 'I': case 'i':
270       fGainFMD3i = MakeSummaryHistogram("gain", "Gains", det, ring);
271       fChi2FMD3i = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
272       break;
273     case 'O': case 'o':
274       fGainFMD3o = MakeSummaryHistogram("gain", "Gains", det, ring);
275       fChi2FMD3o = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
276       break;
277     }
278     break;
279   }
280 }
281
282 //_____________________________________________________________________
283 void AliFMDGainDA::Analyse(UShort_t det, 
284                            Char_t   ring, 
285                            UShort_t sec, 
286                            UShort_t strip) 
287 {
288   // Analyse result of a single strip
289   // 
290   // Parameters: 
291   //  det            Detector number
292   //  ring           Ring identifier 
293   //  sec            Sector number
294   //  strip          Strip number
295   TGraphErrors* grChannel = GetChannel(det,ring,sec,strip);
296   if(!grChannel->GetN()) {
297     AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
298                      det, ring , sec, strip));
299     return;
300   }
301   TF1* fitFunc = new TF1("fitFunc","pol1",-10,280); 
302   fitFunc->SetParameters(100,3);
303   grChannel->Fit(fitFunc,"Q0","",0,fHighPulse);
304   if (grChannel->GetListOfFunctions())
305     grChannel->GetListOfFunctions()->Remove(fitFunc);
306   gROOT->GetListOfFunctions()->Remove(fitFunc);
307
308      
309   Float_t gain    = -1;
310   Float_t error   = -1; 
311   Float_t chi2ndf = -1;
312   if((fitFunc->GetParameter(1)) == (fitFunc->GetParameter(1))) {
313     gain    = fitFunc->GetParameter(1);
314     error   = fitFunc->GetParError(1);
315     if(fitFunc->GetNDF())
316       chi2ndf = fitFunc->GetChisquare() / fitFunc->GetNDF();
317   }
318
319   fOutputFile << det                         << ','
320               << ring                        << ','
321               << sec                         << ','
322               << strip                       << ','
323               << gain                        << ','
324               << error                       << ','
325               << chi2ndf                     <<"\n";
326   
327   //due to RCU trouble, first strips on VAs are excluded
328   // if(strip%128 != 0) {
329     
330   fSummaryGains.SetBinContent(fCurrentSummaryStrip,fitFunc->GetParameter(1));
331   fSummaryGains.SetBinError(fCurrentSummaryStrip,fitFunc->GetParError(1));
332   
333   fCurrentSummaryStrip++;
334
335   TH2* hGain = 0;
336   TH2* hChi2 = 0;
337   switch (det) { 
338   case 1: hGain = fGainFMD1i; hChi2 = fChi2FMD1i; break;
339   case 2: 
340     switch (ring) { 
341     case 'I':  hGain = fGainFMD2i; hChi2 = fChi2FMD2i; break;
342     case 'O':  hGain = fGainFMD2o; hChi2 = fChi2FMD2o; break;
343     }
344     break;
345   case 3:
346     switch (ring) { 
347     case 'I':  hGain = fGainFMD3i; hChi2 = fChi2FMD3i; break;
348     case 'O':  hGain = fGainFMD3o; hChi2 = fChi2FMD3o; break;
349     }
350     break;
351   }
352   if (hGain && hChi2) {
353     Int_t bin = hGain->FindBin(sec, strip);
354     hGain->SetBinContent(bin, gain);
355     hGain->SetBinError(bin, error);
356     hChi2->SetBinContent(bin, chi2ndf);
357   }
358
359     // }
360   if(fSaveHistograms) {
361     TH1F* summary = GetSectorSummary(det, ring, sec);
362     summary->SetBinContent(strip+1, fitFunc->GetParameter(1));
363     summary->SetBinError(strip+1, fitFunc->GetParError(1));
364   }  
365   delete fitFunc;
366 }
367
368 //_____________________________________________________________________
369 void AliFMDGainDA::Terminate(TFile* diagFile)
370 {
371   // End of file 
372   // 
373   // Parameters: 
374   //   None
375   if(diagFile) {
376     diagFile->cd();
377     fSummaryGains.Write();
378   }
379 }
380
381 //_____________________________________________________________________
382 void AliFMDGainDA::WriteHeaderToFile() 
383 {
384   // Write header to the output file
385   // 
386   // Parameters: 
387   //   None
388   AliFMDParameters* pars       = AliFMDParameters::Instance();
389   fOutputFile.write(Form("# %s \n",pars->GetGainShuttleID()),9);
390   TDatime now;
391   fOutputFile << "# This file created from run # " << fRunno 
392               << " @ " << now.AsString() << std::endl;
393   fOutputFile.write("# Detector, "
394                     "Ring, "
395                     "Sector, "
396                     "Strip, "
397                     "Gain, "
398                     "Error, "
399                     "Chi2/NDF \n",56);  
400 }
401
402 //_____________________________________________________________________
403 TH1S* AliFMDGainDA::GetChannelHistogram(UShort_t det, 
404                                         Char_t   ring, 
405                                         UShort_t sec, 
406                                         UShort_t va) 
407 {
408   // Get the current histogram of a single strip
409   // 
410   // Parameters: 
411   //  det            Detector number
412   //  ring           Ring identifier 
413   //  sec            Sector number
414   //  va             VA number
415   Array* secArray  = GetSectorArray(det, ring, sec);
416   Int_t      n         = secArray->GetEntriesFast();
417   Array* cache     = static_cast<Array*>(secArray->At(n-1));
418   TH1S* hChannel       = static_cast<TH1S*>(cache->At(va));
419   
420   return hChannel;
421 }
422
423 //_____________________________________________________________________
424 TGraphErrors* AliFMDGainDA::GetChannel(UShort_t det, 
425                                        Char_t   ring, 
426                                        UShort_t sec, 
427                                        UShort_t strip) 
428 {  
429   // Get the graph of a single strip
430   // 
431   // Parameters: 
432   //  det            Detector number
433   //  ring           Ring identifier 
434   //  sec            Sector number
435   //  strip          Strip number
436   Array*    stripArray = GetStripArray(det, ring, sec, strip);
437   TGraphErrors* hChannel   = static_cast<TGraphErrors*>(stripArray->At(0));
438   
439   return hChannel;
440 }
441
442 //_____________________________________________________________________
443 TH1F* AliFMDGainDA::GetSectorSummary(UShort_t det, 
444                                      Char_t   ring, 
445                                      UShort_t sec) 
446 {
447   Array* secArray    = GetSectorArray(det, ring, sec);
448   Int_t      n           = secArray->GetEntriesFast();
449   return static_cast<TH1F*>(secArray->At(n-2)); // Cache added later
450 }
451
452 //_____________________________________________________________________
453 void AliFMDGainDA::UpdatePulseAndADC(UShort_t det, 
454                                      Char_t ring, 
455                                      UShort_t sec, 
456                                      UShort_t strip) 
457 {
458   // Go to next pulse size
459   // 
460   // Parameters: 
461   //  det            Detector number
462   //  ring           Ring identifier 
463   //  sec            Sector number
464   //  strip          Strip number
465   
466   AliFMDParameters* pars     = AliFMDParameters::Instance();
467   UShort_t          board    = pars->GetAltroMap()->Sector2Board(ring, sec);
468   Int_t             halfring = GetHalfringIndex(det,ring,board/16);
469   
470   if(GetCurrentEvent()> (fNumberOfStripsPerChip*fEventsPerChannel.At(halfring)))
471     return;
472   
473   if ((sec       % 2) && ((strip+1) % fNumberOfStripsPerChip)) return;
474   if (((sec + 1) % 2) && (strip     % fNumberOfStripsPerChip)) return;
475   
476   if(((GetCurrentEvent()) % fPulseLength.At(halfring)) 
477      && GetCurrentEvent() > 0) return;
478      
479   Int_t vaChip   = strip/fNumberOfStripsPerChip; 
480   TH1S* hChannel = GetChannelHistogram(det,ring,sec,vaChip);
481   
482   if(!hChannel->GetEntries()) {
483     AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
484                     det, ring , sec, strip));
485     return;
486   }
487   Double_t mean      = hChannel->GetMean();
488   Double_t rms       = hChannel->GetRMS();
489   Double_t pulse     = (Double_t(fCurrentPulse.At(halfring)) 
490                         * fPulseSize.At(halfring));
491   Int_t    firstBin  = hChannel->GetXaxis()->GetFirst();
492   Int_t    lastBin   = hChannel->GetXaxis()->GetLast();
493   hChannel->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
494   
495   mean               = hChannel->GetMean();
496   rms                = hChannel->GetRMS();
497   
498   hChannel->GetXaxis()->SetRange(firstBin,lastBin);
499   
500   Int_t         channelOff = ((GetCurrentEvent()-1)
501                               / ((fPulseLength.At(halfring)*fHighPulse)
502                                  / fPulseSize.At(halfring)));
503   Int_t         channelNo  = (strip + ((sec % 2 == 1) ? -1 : 1) * channelOff);
504   TGraphErrors* channel    = GetChannel(det,ring,sec,channelNo);
505
506   channel->SetPoint(fCurrentPulse.At(halfring),pulse,mean);
507   channel->SetPointError(fCurrentPulse.At(halfring),0,rms);
508   
509   if(fSaveHistograms) {
510     TH1S* out = 
511       static_cast<TH1S*>(hChannel->Clone(Form("FMD%d%c[%02d,%03d]_0x%02x",
512                                               det, ring, sec, channelNo, 
513                                               int(pulse))));
514     out->SetTitle(Form("FMD%d%c[%02d,%03d] DAC=0x%02x (%3d)",
515                        det, ring, sec, channelNo, int(pulse), int(pulse)));
516     out->SetDirectory(0);
517     Array* arr = GetStripArray(det, ring, sec, channelNo);
518     arr->AddAtAndExpand(out, fCurrentPulse.At(halfring)+1);
519   }
520     
521   hChannel->Reset();
522   
523 }
524
525 //_____________________________________________________________________
526 void AliFMDGainDA::ResetPulseAndUpdateChannel() 
527 {  
528   // Reset all
529   // 
530   // Parameters: 
531   //   None
532   fCurrentPulse.Reset(0); 
533 }
534
535 //_____________________________________________________________________
536 void AliFMDGainDA::FinishEvent() 
537 {
538   // End of event
539   // 
540   // Parameters: 
541   //   None
542   for(UShort_t det=1; det<=3;det++) {
543     UShort_t firstring = (det == 1 ? 1 : 0);
544     for(UShort_t iring = firstring; iring <=1;iring++) {
545       Char_t ring = (iring == 1 ? 'I' : 'O');
546       for(UShort_t board =0 ; board <=1; board++) {
547         Int_t idx = GetHalfringIndex(det,ring,board);
548         
549         if(!fPulseLength.At(idx) || !fEventsPerChannel.At(idx))
550           continue;
551         if(GetCurrentEvent()>0 && 
552            ((GetCurrentEvent() % fPulseLength.At(idx)) == 0))
553           fCurrentPulse.AddAt(fCurrentPulse.At(idx)+1,idx);
554         
555         if(GetCurrentEvent()>0 && 
556            ((GetCurrentEvent()) % fEventsPerChannel.At(idx)) == 0)
557           fCurrentPulse.AddAt(0,idx);
558       }
559     }
560   }
561 }
562 //_____________________________________________________________________
563 //
564 //EOF
565 //