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