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