This is a backward incompatible change in AliRoot. The following methods have been...
[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 "TH1.h"
37 #include "TMath.h"
38 #include "TGraphErrors.h"
39 #include "AliFMDParameters.h"
40 #include "AliFMDAltroMapping.h"
41
42 //_____________________________________________________________________
43 ClassImp(AliFMDGainDA)
44 #if 0 // Do not delete - here to let Emacs indent properly
45 ;
46 #endif
47
48 //_____________________________________________________________________
49 AliFMDGainDA::AliFMDGainDA() 
50   : AliFMDBaseDA(),
51     fGainArray(),
52     fHighPulse(256), 
53     fEventsPerChannel(10),
54     fCurrentPulse(10),
55     fCurrentChannel(10),
56     fNumberOfStripsPerChip(128),
57     fSummaryGains("GainsSummary","Summary of gains",51200,0,51200),
58     fCurrentSummaryStrip(1)
59 {
60   fCurrentPulse.Reset(0);
61   fCurrentChannel.Reset(0);
62   fOutputFile.open("gains.csv");
63   fGainArray.SetOwner(); 
64 }
65
66 //_____________________________________________________________________
67 AliFMDGainDA::AliFMDGainDA(const AliFMDGainDA & gainDA) 
68   :  AliFMDBaseDA(gainDA),
69      fGainArray(gainDA.fGainArray),
70      fHighPulse(gainDA.fHighPulse),
71      fEventsPerChannel(gainDA.fEventsPerChannel),
72      fCurrentPulse(gainDA.fCurrentPulse),
73      fCurrentChannel(gainDA.fCurrentChannel),
74      fNumberOfStripsPerChip(gainDA.fNumberOfStripsPerChip),
75      fSummaryGains(gainDA.fSummaryGains),
76      fCurrentSummaryStrip(gainDA.fCurrentSummaryStrip)
77 {  
78   fCurrentPulse.Reset(0);
79   fCurrentChannel.Reset(0);
80 }
81
82 //_____________________________________________________________________
83 AliFMDGainDA::~AliFMDGainDA() 
84 {
85 }
86
87 //_____________________________________________________________________
88 void AliFMDGainDA::Init() 
89 {
90   
91  
92   Int_t nEventsRequired = 0;
93   
94   //for(UShort_t det=1; det<=3;det++) {
95   //  UShort_t firstring = (det == 1 ? 1 : 0);
96   //  for(UShort_t iring = firstring; iring <=1;iring++) {
97   //    Char_t ring = (iring == 1 ? 'I' : 'O');
98   //    for(UShort_t board =0 ; board <=1; board++) {
99   //    Int_t idx = GetHalfringIndex(det,ring,board);
100   for(Int_t idx = 0;idx<fEventsPerChannel.GetSize();idx++)
101     {
102       
103       Int_t nEvents = 0;
104       if(fPulseSize.At(idx))
105         nEvents = (fPulseLength.At(idx)*fHighPulse) / fPulseSize.At(idx);
106       fEventsPerChannel.AddAt(nEvents,idx);
107       if(nEvents>nEventsRequired) nEventsRequired = nEvents * fNumberOfStripsPerChip;
108       
109     }
110   //}
111   // }
112   
113   //8 pulser values * 128 strips * 100 samples
114   
115   
116   SetRequiredEvents(nEventsRequired); 
117   
118   TObjArray* detArray;
119   TObjArray* ringArray;
120   TObjArray* sectorArray;
121   
122   for(UShort_t det=1;det<=3;det++) {
123     detArray = new TObjArray();
124     detArray->SetOwner();
125     fGainArray.AddAtAndExpand(detArray,det);
126     for (UShort_t ir = 0; ir < 2; ir++) {
127       Char_t   ring = (ir == 0 ? 'O' : 'I');
128       UShort_t nsec = (ir == 0 ? 40  : 20);
129       UShort_t nstr = (ir == 0 ? 2 : 4);
130       ringArray = new TObjArray();
131       ringArray->SetOwner();
132       detArray->AddAtAndExpand(ringArray,ir);
133       for(UShort_t sec =0; sec < nsec;  sec++)  {
134         sectorArray = new TObjArray();
135         sectorArray->SetOwner();
136         ringArray->AddAtAndExpand(sectorArray,sec);
137         for(UShort_t strip = 0; strip < nstr; strip++) {
138           TH1S* hChannel = new TH1S(Form("hFMD%d%c_%d_%d",det,ring,sec,strip),
139                                     Form("hFMD%d%c_%d_%d",det,ring,sec,strip),
140                                     1024,0,1023);
141           hChannel->SetDirectory(0);
142           sectorArray->AddAtAndExpand(hChannel,strip);
143         }
144       }
145     }
146   }
147 }
148
149 //_____________________________________________________________________
150 void AliFMDGainDA::AddChannelContainer(TObjArray* sectorArray, 
151                                        UShort_t det  , 
152                                        Char_t   ring,  
153                                        UShort_t sec, 
154                                        UShort_t strip) 
155 {  
156   TGraphErrors* hChannel  = new TGraphErrors();
157   hChannel->SetName(Form("FMD%d%c[%02d,%03d]", det, ring, sec, strip));
158   hChannel->SetTitle(Form("FMD%d%c[%02d,%03d] ADC vs DAC", 
159                           det, ring, sec, strip));
160   sectorArray->AddAtAndExpand(hChannel,strip);
161 }
162
163 //_____________________________________________________________________
164 void AliFMDGainDA::FillChannels(AliFMDDigit* digit) {
165
166   UShort_t det   = digit->Detector();
167   Char_t   ring  = digit->Ring();
168   UShort_t sec   = digit->Sector();
169   UShort_t strip = digit->Strip();
170   
171   //Strip is always seen as the first in a VA chip. All other strips are junk.
172   //Strips are counted from zero on even sectors and from 511 on odd sectors...
173    
174   if((sec%2)     && ((strip+1) % fNumberOfStripsPerChip)) return;
175   if(((sec+1)%2) && (strip % fNumberOfStripsPerChip)) return;
176   
177   Int_t vaChip   = strip / fNumberOfStripsPerChip; 
178   TH1S* hChannel = GetChannelHistogram(det, ring, sec, vaChip);
179   hChannel->Fill(digit->Counts());
180   UpdatePulseAndADC(det,ring,sec,strip);
181 }
182
183 //_____________________________________________________________________
184 void AliFMDGainDA::Analyse(UShort_t det, 
185                            Char_t   ring, 
186                            UShort_t sec, 
187                            UShort_t strip) {
188   TGraphErrors* grChannel = GetChannel(det,ring,sec,strip);
189   if(!grChannel->GetN()) {
190     AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
191                      det, ring , sec, strip));
192     return;
193   }
194   TF1 fitFunc("fitFunc","pol1",-10,280); 
195   fitFunc.SetParameters(100,3);
196   grChannel->Fit("fitFunc","Q0+","",0,fHighPulse);
197      
198   Float_t gain    = -1;
199   Float_t error   = -1; 
200   Float_t chi2ndf = -1;
201   if((fitFunc.GetParameter(1)) == (fitFunc.GetParameter(1))) {
202     gain    = fitFunc.GetParameter(1);
203     error   = fitFunc.GetParError(1);
204     if(fitFunc.GetNDF())
205       chi2ndf = fitFunc.GetChisquare() / fitFunc.GetNDF();
206   }
207   
208   fOutputFile << det                         << ','
209               << ring                        << ','
210               << sec                         << ','
211               << strip                       << ','
212               << gain                        << ','
213               << error                       << ','
214               << chi2ndf                     <<"\n";
215   
216   //due to RCU trouble, first strips on VAs are excluded
217   // if(strip%128 != 0) {
218     
219   fSummaryGains.SetBinContent(fCurrentSummaryStrip,fitFunc.GetParameter(1));
220   fSummaryGains.SetBinError(fCurrentSummaryStrip,fitFunc.GetParError(1));
221   
222   fCurrentSummaryStrip++;
223     // }
224   if(fSaveHistograms) {
225     gDirectory->cd(GetSectorPath(det,ring, sec, kTRUE));
226     
227     TH1F* summary = dynamic_cast<TH1F*>(gDirectory->Get("Summary"));
228     if (!summary) { 
229       Int_t nStr = (ring == 'I' ? 512 : 256);
230       summary = new TH1F("Summary", Form("Summary of gains in FMD%d%c[%02d]", 
231                                          det, ring, sec), 
232                          nStr, -.5, nStr-.5);
233       summary->SetXTitle("Strip");
234       summary->SetYTitle("Gain [ADC/DAC]");
235       summary->SetDirectory(gDirectory);
236     }
237     summary->SetBinContent(strip+1, fitFunc.GetParameter(1));
238     summary->SetBinError(strip+1, fitFunc.GetParError(1));
239     
240     gDirectory->cd(GetStripPath(det,ring,sec,strip, kTRUE));
241     grChannel->SetName(Form("FMD%d%c[%02d,%03d]",det,ring,sec,strip));
242     // grChannel->SetDirectory(gDirectory);
243     grChannel->Write();
244     // grChannel->Write(Form("grFMD%d%c_%d_%d",det,ring,sec,strip));
245   }  
246 }
247
248 //_____________________________________________________________________
249 void AliFMDGainDA::Terminate(TFile* diagFile)
250 {
251   if(diagFile) {
252     diagFile->cd();
253     fSummaryGains.Write();
254   }
255 }
256
257 //_____________________________________________________________________
258 void AliFMDGainDA::WriteHeaderToFile() 
259 {
260   AliFMDParameters* pars       = AliFMDParameters::Instance();
261   fOutputFile.write(Form("# %s \n",pars->GetGainShuttleID()),9);
262   fOutputFile.write("# Detector, "
263                     "Ring, "
264                     "Sector, "
265                     "Strip, "
266                     "Gain, "
267                     "Error, "
268                     "Chi2/NDF \n",56);
269   
270 }
271
272 //_____________________________________________________________________
273 TH1S* AliFMDGainDA::GetChannelHistogram(UShort_t det, 
274                                         Char_t   ring, 
275                                         UShort_t sec, 
276                                         UShort_t strip) 
277 {
278   
279   UShort_t  Ring = 1;
280   if(ring == 'O')
281     Ring = 0;
282   
283   
284   TObjArray* detArray  = static_cast<TObjArray*>(fGainArray.At(det));
285   TObjArray* ringArray = static_cast<TObjArray*>(detArray->At(Ring));
286   TObjArray* secArray  = static_cast<TObjArray*>(ringArray->At(sec));
287   TH1S* hChannel       = static_cast<TH1S*>(secArray->At(strip));
288   
289   return hChannel;
290 }
291
292 //_____________________________________________________________________
293 TGraphErrors* AliFMDGainDA::GetChannel(UShort_t det, 
294                                        Char_t   ring, 
295                                        UShort_t sec, 
296                                        UShort_t strip) 
297 {  
298   UShort_t      iring     = (ring == 'O' ? 0 : 1);
299   TObjArray*    detArray  = static_cast<TObjArray*>(fDetectorArray.At(det));
300   TObjArray*    ringArray = static_cast<TObjArray*>(detArray->At(iring));
301   TObjArray*    secArray  = static_cast<TObjArray*>(ringArray->At(sec));
302   TGraphErrors* hChannel  = static_cast<TGraphErrors*>(secArray->At(strip));
303   
304   return hChannel;
305 }
306
307 //_____________________________________________________________________
308 void AliFMDGainDA::UpdatePulseAndADC(UShort_t det, 
309                                      Char_t ring, 
310                                      UShort_t sec, 
311                                      UShort_t strip) 
312 {
313   
314   AliFMDParameters* pars = AliFMDParameters::Instance();
315   // UInt_t ddl, board,chip,ch;
316   UShort_t board = pars->GetAltroMap()->Sector2Board(ring, sec);
317   // pars->Detector2Hardware(det,ring,sec,strip,ddl,board,chip,ch);
318   /// pars->GetAltroMap()->Strip2Channel(
319   Int_t halfring = GetHalfringIndex(det,ring,board/16);
320   
321   if(GetCurrentEvent()> (fNumberOfStripsPerChip*fEventsPerChannel.At(halfring)))
322     return;
323   
324   if((sec%2)     && ((strip+1) % fNumberOfStripsPerChip)) return;
325   
326   if(((sec+1)%2) && (strip % fNumberOfStripsPerChip)) return;
327   
328   if(((GetCurrentEvent()) % fPulseLength.At(halfring)) 
329      && GetCurrentEvent() > 0) return;
330      
331   Int_t vaChip = strip/fNumberOfStripsPerChip; 
332   TH1S* hChannel = GetChannelHistogram(det,ring,sec,vaChip);
333   
334   if(!hChannel->GetEntries()) {
335     AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
336                     det, ring , sec, strip));
337     return;
338   }
339   Double_t mean      = hChannel->GetMean();
340   Double_t rms       = hChannel->GetRMS();
341   Double_t pulse     = (Double_t(fCurrentPulse.At(halfring)) 
342                         * fPulseSize.At(halfring));
343   Int_t    firstBin  = hChannel->GetXaxis()->GetFirst();
344   Int_t    lastBin   = hChannel->GetXaxis()->GetLast();
345   hChannel->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
346   
347   mean               = hChannel->GetMean();
348   rms                = hChannel->GetRMS();
349   
350   hChannel->GetXaxis()->SetRange(firstBin,lastBin);
351   
352   Int_t    channelNumber      = (strip + 
353                                  (GetCurrentEvent()-1)
354                                  / ((fPulseLength.At(halfring)*fHighPulse)
355                                     / fPulseSize.At(halfring))); 
356   if(sec%2)
357     channelNumber      = (strip - 
358                           (GetCurrentEvent()-1)
359                           / ((fPulseLength.At(halfring)*fHighPulse)
360                              / fPulseSize.At(halfring))); 
361   
362   TGraphErrors* channel = GetChannel(det,ring,sec,channelNumber);
363   
364   channel->SetPoint(fCurrentPulse.At(halfring),pulse,mean);
365   channel->SetPointError(fCurrentPulse.At(halfring),0,rms);
366   
367   if(fSaveHistograms) {
368     gDirectory->cd(GetStripPath(det,ring,sec,channelNumber));
369     hChannel->Write(Form("%s_pulse_%03d",hChannel->GetName(),(Int_t)pulse));
370     
371   }
372     
373   hChannel->Reset();
374   
375 }
376
377 //_____________________________________________________________________
378 void AliFMDGainDA::ResetPulseAndUpdateChannel() 
379 {  
380   fCurrentPulse.Reset(0); 
381 }
382
383 //_____________________________________________________________________
384 void AliFMDGainDA::FinishEvent() 
385 {
386   for(UShort_t det=1; det<=3;det++) {
387     UShort_t firstring = (det == 1 ? 1 : 0);
388     for(UShort_t iring = firstring; iring <=1;iring++) {
389       Char_t ring = (iring == 1 ? 'I' : 'O');
390       for(UShort_t board =0 ; board <=1; board++) {
391         Int_t idx = GetHalfringIndex(det,ring,board);
392         
393         if( !fPulseLength.At(idx) || !fEventsPerChannel.At(idx))
394           continue;
395         if(GetCurrentEvent()>0 && ((GetCurrentEvent() % fPulseLength.At(idx)) == 0))
396           fCurrentPulse.AddAt(fCurrentPulse.At(idx)+1,idx);
397         
398         if(GetCurrentEvent()>0 && ((GetCurrentEvent()) % fEventsPerChannel.At(idx)) == 0)
399           fCurrentPulse.AddAt(0,idx);
400       }
401     }
402   }
403 }
404 //_____________________________________________________________________
405 //
406 //EOF
407 //