Removed a std::cout statement
[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 // This class contains the implementation of the gain detector algorithms (DA) for the FMD.
22 // The data is collected in histograms that are reset for each pulse length after the 
23 // mean and standard deviation are put into a TGraphErrors object. After a certain number of pulses
24 // (usually 8) the graph is fitted to a straight line. The gain is then slope of this line as
25 // it combines the known pulse and the response of the detector.
26
27
28 #include "AliFMDGainDA.h"
29 #include "iostream"
30 #include "fstream"
31 #include "AliLog.h"
32 #include "TF1.h"
33 #include "TH1.h"
34 #include "TMath.h"
35 #include "TGraphErrors.h"
36 #include "AliFMDParameters.h"
37
38 //_____________________________________________________________________
39 ClassImp(AliFMDGainDA)
40
41 //_____________________________________________________________________
42 AliFMDGainDA::AliFMDGainDA() : AliFMDBaseDA(),
43   fGainArray(),
44   fPulseSize(32),
45   fHighPulse(256), 
46   fPulseLength(100),
47   fEventsPerChannel(0),
48   fCurrentPulse(0),
49   fCurrentChannel(0),
50   fNumberOfStripsPerChip(128)
51 {
52   fOutputFile.open("gains.csv");
53   fGainArray.SetOwner();
54   
55 }
56
57 //_____________________________________________________________________
58 AliFMDGainDA::AliFMDGainDA(const AliFMDGainDA & gainDA) :  
59   AliFMDBaseDA(gainDA),
60   fGainArray(gainDA.fGainArray),
61   fPulseSize(gainDA.fPulseSize),
62   fHighPulse(gainDA.fHighPulse),
63   fPulseLength(gainDA.fPulseLength),
64   fEventsPerChannel(gainDA.fEventsPerChannel),
65   fCurrentPulse(gainDA.fCurrentPulse),
66   fCurrentChannel(gainDA.fCurrentChannel),
67   fNumberOfStripsPerChip(gainDA.fNumberOfStripsPerChip)
68 {
69   
70 }
71
72 //_____________________________________________________________________
73 AliFMDGainDA::~AliFMDGainDA() {
74   
75  
76
77 }
78
79 //_____________________________________________________________________
80 void AliFMDGainDA::Init() {
81     
82   fEventsPerChannel = (fPulseLength*fHighPulse) / fPulseSize ;
83   
84   SetRequiredEvents(fEventsPerChannel*fNumberOfStripsPerChip); //8 pulser values * 128 strips * 100 samples
85   TObjArray* detArray;
86   TObjArray* ringArray;
87   TObjArray* sectorArray;
88   
89   for(UShort_t det=1;det<=3;det++) {
90     detArray = new TObjArray();
91     detArray->SetOwner();
92     fGainArray.AddAtAndExpand(detArray,det);
93     for (UShort_t ir = 0; ir < 2; ir++) {
94       Char_t   ring = (ir == 0 ? 'O' : 'I');
95       UShort_t nsec = (ir == 0 ? 40  : 20);
96       UShort_t nstr = (ir == 0 ? 2 : 4);
97       ringArray = new TObjArray();
98       ringArray->SetOwner();
99       detArray->AddAtAndExpand(ringArray,ir);
100       for(UShort_t sec =0; sec < nsec;  sec++)  {
101         sectorArray = new TObjArray();
102         sectorArray->SetOwner();
103         ringArray->AddAtAndExpand(sectorArray,sec);
104         for(UShort_t strip = 0; strip < nstr; strip++) {
105           TH1S* hChannel = new TH1S(Form("hFMD%d%c_%d_%d",det,ring,sec,strip),Form("hFMD%d%c_%d_%d",det,ring,sec,strip),1024,0,1023);
106           hChannel->SetDirectory(0);
107           sectorArray->AddAtAndExpand(hChannel,strip);
108         }
109       }
110     }
111   }
112 }
113
114 //_____________________________________________________________________
115 void AliFMDGainDA::AddChannelContainer(TObjArray* sectorArray, 
116                                        UShort_t det, 
117                                        Char_t ring, 
118                                        UShort_t sec, 
119                                        UShort_t strip) {
120   
121   TGraphErrors* hChannel  = new TGraphErrors();
122   sectorArray->AddAtAndExpand(hChannel,strip);
123 }
124
125 //_____________________________________________________________________
126 void AliFMDGainDA::FillChannels(AliFMDDigit* digit) {
127
128   UShort_t det   = digit->Detector();
129   Char_t   ring  = digit->Ring();
130   UShort_t sec   = digit->Sector();
131   UShort_t strip = digit->Strip();
132   
133   
134   if(strip%fNumberOfStripsPerChip) return;
135   Int_t VAchip = strip / fNumberOfStripsPerChip; 
136   TH1S* hChannel = GetChannelHistogram(det,ring,sec,VAchip);
137   hChannel->Fill(digit->Counts());
138   UpdatePulseAndADC(det,ring,sec,strip);
139 }
140
141 //_____________________________________________________________________
142 void AliFMDGainDA::Analyse(UShort_t det, 
143                            Char_t ring, 
144                            UShort_t sec, 
145                            UShort_t strip) {
146   TGraphErrors* grChannel = GetChannel(det,ring,sec,strip);
147   if(!grChannel->GetN()) {
148     AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",det, ring , sec, strip));
149     return;
150   }
151   TF1* fitFunc = new TF1("fitFunc","pol1",-10,280); 
152   fitFunc->SetParameters(100,3);
153   grChannel->Fit("fitFunc","Q","Q",0,fHighPulse);
154   AliFMDParameters* pars = AliFMDParameters::Instance();
155   UInt_t ddl, board,chip,channel;
156   pars->Detector2Hardware(det,ring,sec,strip,ddl,board,chip,channel);
157     
158   Float_t chi2ndf = 0;
159   if(fitFunc->GetNDF())
160     chi2ndf = fitFunc->GetChisquare() / fitFunc->GetNDF();
161   ddl = ddl + kBaseDDL;
162   
163   Int_t relStrip = strip%fNumberOfStripsPerChip;
164   
165   fOutputFile << ddl                         << ','
166               << board                       << ','
167               << chip                        << ','
168               << channel                     << ','
169               << relStrip                    << ','
170               << fitFunc->GetParameter(1)    << ','
171               << fitFunc->GetParError(1)     << ','
172               << chi2ndf                     <<"\n";
173   
174   
175   if(fSaveHistograms) {
176     
177     gDirectory->cd(Form("%s:FMD%d%c/sector_%d/strip_%d",fDiagnosticsFilename,det,ring,sec,strip));
178     grChannel->Write(Form("grFMD%d%c_%d_%d",det,ring,sec,strip));
179   }
180   
181   delete fitFunc;
182   delete grChannel;
183 }
184
185 //_____________________________________________________________________
186 void AliFMDGainDA::WriteHeaderToFile() {
187   fOutputFile.write("# Gains \n",9);
188   fOutputFile.write("# Rcu, Board, Chip, Channel, Strip, Gain, Error, Chi2/NDF \n",59);
189   
190 }
191
192 //_____________________________________________________________________
193 TH1S* AliFMDGainDA::GetChannelHistogram(UShort_t det, 
194                                         Char_t ring, 
195                                         UShort_t sec, 
196                                         UShort_t strip) {
197   
198   UShort_t  Ring = 1;
199   if(ring == 'O')
200     Ring = 0;
201   
202   
203   TObjArray* detArray  = static_cast<TObjArray*>(fGainArray.At(det));
204   TObjArray* ringArray = static_cast<TObjArray*>(detArray->At(Ring));
205   TObjArray* secArray  = static_cast<TObjArray*>(ringArray->At(sec));
206   TH1S* hChannel       = static_cast<TH1S*>(secArray->At(strip));
207   
208   return hChannel;
209 }
210
211 //_____________________________________________________________________
212 TGraphErrors* AliFMDGainDA::GetChannel(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip) {
213   
214   UShort_t  Ring = 1;
215   if(ring == 'O')
216     Ring = 0;
217   
218   
219   TObjArray* detArray          = static_cast<TObjArray*>(fDetectorArray.At(det));
220   TObjArray* ringArray         = static_cast<TObjArray*>(detArray->At(Ring));
221   TObjArray* secArray          = static_cast<TObjArray*>(ringArray->At(sec));
222   TGraphErrors* hChannel       = static_cast<TGraphErrors*>(secArray->At(strip));
223   
224   return hChannel;
225 }
226
227 //_____________________________________________________________________
228 void AliFMDGainDA::UpdatePulseAndADC(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip) {
229   
230   if(strip%fNumberOfStripsPerChip) return;
231   if(((GetCurrentEvent())%fPulseLength) && GetCurrentEvent()>0) return;
232     
233   Int_t VAchip = strip/fNumberOfStripsPerChip; 
234   TH1S* hChannel = GetChannelHistogram(det,ring,sec,VAchip);
235   
236   if(!hChannel->GetEntries()) {
237     AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",det, ring , sec, strip));
238     return;
239   }
240   Double_t mean      = hChannel->GetMean();
241   Double_t rms       = hChannel->GetRMS();
242   Double_t pulse     = (Double_t)fCurrentPulse*fPulseSize;
243   
244   for(Int_t i=1;i<hChannel->GetNbinsX();i++)
245     if(hChannel->GetBinCenter(i)>mean+5*rms || hChannel->GetBinCenter(i)<mean-5*rms)
246       hChannel->SetBinContent(i,0);
247   
248   mean      = hChannel->GetMean();
249   rms       = hChannel->GetRMS();
250     
251   Int_t channelNumber = strip + (GetCurrentEvent()-1)/800; 
252   
253   TGraphErrors* channel = GetChannel(det,ring,sec,channelNumber);
254   
255   channel->SetPoint(fCurrentPulse,pulse,mean);
256   channel->SetPointError(fCurrentPulse,0,rms);
257   
258   if(fSaveHistograms) {
259     gDirectory->cd(Form("%s:FMD%d%c/sector_%d/strip_%d",fDiagnosticsFilename,det,ring,sec,channelNumber));
260     hChannel->Write(Form("hFMD%d%c_%d_%d_pulse_%d",det,ring,sec,channelNumber,fCurrentPulse));
261   }
262   
263   hChannel->Reset();
264   
265 }
266
267 //_____________________________________________________________________
268 void AliFMDGainDA::ResetPulseAndUpdateChannel() {
269   
270   fCurrentPulse   = 0;
271   
272 }
273
274 //_____________________________________________________________________
275 void AliFMDGainDA::FinishEvent() {
276   
277   if(GetCurrentEvent()>0 && (GetCurrentEvent()%fPulseLength==0))
278     fCurrentPulse++;
279   
280   if(GetCurrentEvent()>0 && (GetCurrentEvent())%fEventsPerChannel==0)
281     fCurrentPulse = 0;
282
283
284 }
285 //_____________________________________________________________________
286 //
287 //EOF
288 //