Fixes for coverity
[u/mrichter/AliRoot.git] / PMD / AliPMDCalibrator.cxx
1 //////////////////////////////////////////////////////////////////////////////
2 //
3 // * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 // *                                                                        *
5 // * Author: The ALICE Off-line Project.                                    *
6 // * Contributors are mentioned in the code where appropriate.              *
7 // *                                                                        *
8 // * Permission to use, copy, modify and distribute this software and its   *
9 // * documentation strictly for non-commercial purposes is hereby granted   *
10 // * without fee, provided that the above copyright notice appears in all   *
11 // * copies and that both the copyright notice and this permission notice   *
12 // * appear in the supporting documentation. The authors make no claims     *
13 // * about the suitability of this software for any purpose. It is          *
14 // * provided "as is" without express or implied warranty.                  *
15 // **************************************************************************/
16 //
17 //////////////////////////////////////////////////////////////////////////////
18
19 #include "TF1.h"
20 #include "TFile.h"
21 #include "TObjString.h"
22 #include "TROOT.h"
23 #include "TClonesArray.h"
24 #include "TH1F.h"
25 #include "TObjArray.h"
26 #include "TMath.h"
27
28 // --- Standard library ---
29
30 // --- AliRoot header files ---
31 #include "AliLog.h"
32 #include "AliRawReaderFile.h"
33 #include "AliPMDCalibrator.h"
34 #include "AliRawReaderDate.h"
35 #include "AliPMDRawStream.h"
36 #include "AliPMDCalibData.h"
37 #include "AliPMDddldata.h"
38 #include "AliCDBManager.h"
39 #include "AliCDBId.h"
40 #include "AliCDBMetaData.h"
41 #include "AliDAQ.h"
42
43 //#include "AliPMDCleanNoise.h"
44 //#include "AliPMDCleaner.h"
45
46 #include "AliPMDPedestal.h"
47 #include "AliCDBManager.h"
48 #include "AliCDBEntry.h"
49
50 ClassImp(AliPMDCalibrator)
51
52
53 AliPMDCalibrator::AliPMDCalibrator():
54   fCalibGain(new AliPMDCalibData()), fCalibPed(new AliPMDPedestal())
55 {
56   // Standard Constructor
57   for(Int_t idet = 0; idet < kDet; idet++)
58     {
59       fHdetIso[idet] = NULL ;
60       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
61         {
62           fHsmIso[idet][ismn] = NULL ;
63           for(Int_t jrow = 0; jrow < kMaxRow; jrow++)
64             {
65               for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
66                 {
67                   fGainFact[idet][ismn][jrow][kcol] = 0.0;
68                   fHadcIso[idet][ismn][jrow][kcol]  = NULL;
69                 }
70             }
71         }
72     }
73 }
74 // ------------------------------------------------------------------------ //
75 AliPMDCalibrator::AliPMDCalibrator(const AliPMDCalibrator &pmdcalibrator):
76   fCalibGain(new AliPMDCalibData()), fCalibPed(new AliPMDPedestal())
77 {
78   for(Int_t idet = 0; idet < 2; idet++)
79     {
80       fHdetIso[idet] = pmdcalibrator.fHdetIso[idet] ;
81       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
82         {
83           fHsmIso[idet][ismn] = pmdcalibrator.fHsmIso[idet][ismn] ;
84           for(Int_t jrow = 0; jrow < kMaxRow; jrow++)
85             {
86               for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
87                 {
88                   fGainFact[idet][ismn][jrow][kcol] = pmdcalibrator.fGainFact[idet][ismn][jrow][kcol];
89                   fHadcIso[idet][ismn][jrow][kcol]  = pmdcalibrator.fHadcIso[idet][ismn][jrow][kcol];
90                 }
91             }
92         }
93     }
94   
95 }
96 // ------------------------------------------------------------------------ //
97 AliPMDCalibrator &AliPMDCalibrator::operator=(const AliPMDCalibrator &pmdcalibrator)
98 {
99   if(this != &pmdcalibrator)
100     {
101       for(Int_t idet = 0; idet < kDet; idet++)
102         {
103           fHdetIso[idet] = pmdcalibrator.fHdetIso[idet] ;
104           for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
105             {
106               fHsmIso[idet][ismn] = pmdcalibrator.fHsmIso[idet][ismn] ;
107               for(Int_t jrow = 0; jrow < kMaxRow;jrow++)
108                 {
109                   for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
110                     {
111                       fGainFact[idet][ismn][jrow][kcol] =
112                         pmdcalibrator.fGainFact[idet][ismn][jrow][kcol];
113                       fHadcIso[idet][ismn][jrow][kcol]  =
114                         pmdcalibrator.fHadcIso[idet][ismn][jrow][kcol];
115                     }
116                 }
117             }
118         }
119     }
120   return *this;
121 }
122 // ------------------------------------------------------------------------ //
123 AliPMDCalibrator::~AliPMDCalibrator()
124 {
125   // destructor
126   for (Int_t i=0; i<kDet; i++)
127     delete fHdetIso[i] ;
128   for (Int_t i=0; i<kDet; i++)
129     for (Int_t j=0; j<kMaxSMN; j++)
130       delete fHsmIso[i][j] ;
131   for (Int_t i=0; i<kDet; i++)
132     for (Int_t j=0; j<kMaxSMN; j++)
133       for (Int_t k=0; k<kMaxRow; k++)
134         for (Int_t l=0; l<kMaxCol; l++)
135           delete fHadcIso[i][j][k][l] ;
136   delete fCalibGain;
137   delete fCalibPed;
138 }
139 // ------------------------------------------------------------------------ //
140
141 void AliPMDCalibrator::Exec()
142 {
143   // reads parameters and does the calibration
144   CalculateIsoCell() ;
145 }
146 // ------------------------------------------------------------------------ //
147
148 void AliPMDCalibrator::Init()
149 {
150   // intializes everything
151   char hname[50];
152   char hname2[50];
153   char hname24[50];
154   char hnameiso[120];
155   char htitle1[120];
156   char htitle2[120];
157   
158   for(Int_t d = 0; d < kDet; d++) {
159     snprintf(hname2,50,"Isolated cell adc for Det Plane %d",d);
160     fHdetIso[d]= new TH1F(hname2,htitle2,100,0,1000);
161     for(Int_t i1 = 0; i1 < kMaxSMN; i1++) {
162       snprintf(hname,50,"det_%d_iso_sm_%2d",d,i1);
163       snprintf(hname24,50,"det_%d_iso_sm_%2d",d,i1);
164       fHsmIso[d][i1]= new TH1F(hname,hname24,100,0,1000);
165       for(Int_t j1 = 0; j1 < kMaxRow; j1++) {
166         for(Int_t k1 = 0; k1 < kMaxCol; k1++) {
167           snprintf(hnameiso,120,"Isolated Cell ADC for det_%d_cell_sm%d_row%d_col%d"
168                   ,d,i1,j1,k1);
169           snprintf(htitle1,120,"Isolated Cell ADC for det_%d_cell_sm%d_row%d_col%d"
170                   ,d,i1,j1,k1);
171           
172           TObject *old=gDirectory->GetList()->FindObject(hnameiso);
173           if (old) gDirectory->GetList()->Remove(old);
174           fHadcIso[d][i1][j1][k1] = new TH1F(hnameiso,htitle1,100,0.,1000.);
175         }
176       }
177     }
178   }
179   
180 }
181
182 // ------------------------------------------------------------------------ //
183
184 void AliPMDCalibrator::CalculateIsoCell()
185 {
186   // Calculates the ADC of isolated cell
187
188   TObjArray pmdddlcont;
189
190   const Int_t kCellNeighbour = 6;
191
192   Int_t neibx[6] = {1,0,-1,-1,0,1};
193   Int_t neiby[6] = {0,1,1,0,-1,-1};
194
195   Int_t id1 = 0,jd1 = 0;        //neighbour row/col
196   Int_t countisocell = 0 ;      //number of isilated cell
197   Int_t isocount = 0;           //number of neighbours with 0 signal
198   Int_t maxhit = 0;
199   Int_t nhit[kDet][kMaxSMN];
200   Int_t d1[kDet][kMaxSMN][kMaxRow][kMaxCol];
201   Int_t nhitcell[kDet][kMaxSMN][kMaxRow][kMaxCol];
202   
203   for(Int_t idet = 0; idet < kDet; idet++)
204     {
205       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
206         {
207           nhit[idet][ismn] = 0;
208           for(Int_t irow = 0; irow < kMaxRow; irow++)
209             {
210               for(Int_t icol = 0; icol < kMaxCol; icol++)
211                 {
212                   d1[idet][ismn][irow][icol] = 0;
213                   nhitcell[idet][ismn][irow][icol] = 0;
214                 }
215             }
216         }
217     }
218
219   Float_t tempnhit1[kDet][kMaxSMN];
220   Float_t tempnhit2[kDet][kMaxSMN];
221   Float_t meannhit[kDet][kMaxSMN];
222   Float_t meanSqnhit[kDet][kMaxSMN];
223   Float_t sigmanhit[kDet][kMaxSMN];
224   Float_t count[kDet][kMaxSMN];
225   Float_t nhitcut[kDet][kMaxSMN];
226
227   for (Int_t idet = 0; idet < kDet; idet++)
228   {
229       for (Int_t ismn = 0; ismn < kMaxSMN; ismn++)
230       {
231           tempnhit1[idet][ismn]  = 0.;
232           tempnhit2[idet][ismn]  = 0.;
233           meannhit[idet][ismn]   = 0.;
234           meanSqnhit[idet][ismn] = 0.;
235           sigmanhit[idet][ismn]  = 0.;
236           count[idet][ismn]      = 0.;
237           nhitcut[idet][ismn]    = 0.;
238       }
239   }
240
241
242   //accessing raw data
243   AliRawReaderFile reader(".");
244   AliPMDRawStream stream(&reader);
245   while(reader.NextEvent())
246     { 
247       // New PMD Reader is plugged in
248         Int_t iddl = -1;
249         while ((iddl = stream.DdlData(&pmdddlcont)) >=0) {
250
251             //reader.Select("PMD", iddl, iddl);
252             //stream.DdlData(iddl,&pmdddlcont);
253           Int_t ientries = pmdddlcont.GetEntries();
254           for (Int_t ient = 0; ient < ientries; ient++)
255             {
256               AliPMDddldata *pmdddl = 
257                 (AliPMDddldata*)pmdddlcont.UncheckedAt(ient);
258               Int_t idet = pmdddl->GetDetector();
259               Int_t ismn = pmdddl->GetSMN();
260               Int_t irow = pmdddl->GetRow();
261               Int_t icol = pmdddl->GetColumn();
262               Float_t isig1 = pmdddl->GetSignal();
263               // Pedestal Subtraction
264               Int_t   pedmeanrms = 
265                 fCalibPed->GetPedMeanRms(idet,ismn,irow,icol);
266               Int_t   pedrms1    = (Int_t) pedmeanrms%100;
267               Float_t pedrms     = (Float_t)pedrms1/10.;
268               Float_t pedmean    = (Float_t) (pedmeanrms - pedrms1)/1000.0;
269               Float_t isig = isig1 - (pedmean + 3.0*pedrms);
270               if (isig>0)
271                 {
272                   d1[idet][ismn][irow][icol] = (Int_t)isig;
273                   nhitcell[idet][ismn][irow][icol] += 1; 
274                 }
275             }//ient loop
276           pmdddlcont.Clear();
277         }//iddl loop
278       maxhit = 0;
279       for(Int_t idet=0; idet < kDet; idet++)
280         {
281           for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
282             {
283               for(Int_t irow = 0; irow < kMaxRow; irow++)
284                 {
285                   for(Int_t icol = 0; icol < kMaxCol; icol++)
286                     {
287                       if(d1[idet][ismn][irow][icol] > 0)
288                         {
289                           isocount = 0;
290                           for(Int_t ii = 0; ii < kCellNeighbour; ii++)
291                             {
292                               id1 = irow + neibx[ii];
293                               jd1 = icol + neiby[ii];
294                               if(d1[idet][ismn][id1][jd1] == 0)
295                                 {
296                                   isocount++;
297                                   if(isocount == kCellNeighbour)
298                                     {
299                                       countisocell++;
300                                       maxhit++;
301                                       fHdetIso[idet]->
302                                         Fill(d1[idet][ismn][irow][icol]);
303                                       fHsmIso[idet][ismn]->
304                                         Fill(d1[idet][ismn][irow][icol]);
305                                       fHadcIso[idet][ismn][irow][icol]->
306                                         Fill(d1[idet][ismn][irow][icol]);
307                                     }
308                                 }
309                             }  // neigh cell cond.
310                         }     // d>0 cond.
311                     }
312                 }
313             }
314         } //det
315     }//event loop
316
317   // Mean and Sigma Calculations
318   for(Int_t idet=0; idet < kDet; idet++){
319       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++){
320           for(Int_t irow = 0; irow < kMaxRow; irow++){
321               for(Int_t icol = 0; icol < kMaxCol; icol++){
322                   if(nhitcell[idet][ismn][irow][icol]>0){
323                       count[idet][ismn] += 1;
324                       tempnhit1[idet][ismn] += nhitcell[idet][ismn][irow][icol];
325                       tempnhit2[idet][ismn] += nhitcell[idet][ismn][irow][icol]
326                           *nhitcell[idet][ismn][irow][icol];
327                   }
328               }
329           }
330       }
331   }//det loop
332
333   //cout<<"nhit cell = "<<idet<<"  "<<ismn<<"  "
334   //  <<irow<<"  "<<icol<<"  "<<nhitcell[idet][ismn][irow][icol]<<endl;
335   //count[idet][ismn] += 1;
336
337   for(Int_t i=0; i < kDet; i++)
338     {
339       for(Int_t j=0; j < kMaxSMN; j++)
340         {
341           if(count[i][j] > 0.)
342             {
343                 meannhit[i][j]   = tempnhit1[i][j]/count[i][j];
344                 meanSqnhit[i][j] = tempnhit2[i][j]/count[i][j]; 
345                 sigmanhit[i][j]  = sqrt(meanSqnhit[i][j]-
346                                       (meannhit[i][j]*meannhit[i][j]));
347                 nhitcut[i][j]    = 3*sigmanhit[i][j] + meannhit[i][j];
348             }
349         }
350     }
351
352
353   Double_t histdetMean[kDet];
354   Double_t histMean[kDet][kMaxSMN];
355   Double_t isoMean[kDet][kMaxSMN][kMaxRow][kMaxCol];
356   Double_t smNormFactor[kDet][kMaxSMN];
357   
358   for(Int_t det1 = 0; det1 < kDet; det1++)
359     {
360       histdetMean[det1]= fHdetIso[det1]->GetMean();
361       for(Int_t i1 = 0; i1 < kMaxSMN; i1++)
362         {
363           histMean[det1][i1]= fHsmIso[det1][i1]->GetMean();
364           if(histMean[det1][i1]>0.0 && histdetMean[det1]>0.0)
365             {
366               smNormFactor[det1][i1]= histdetMean[det1]/histMean[det1][i1];
367             }
368           for(Int_t j1 = 0; j1 < kMaxRow; j1++)
369             {
370               for(Int_t k1 = 0; k1 < kMaxCol; k1++)
371                 {
372                  if(nhitcell[det1][i1][j1][k1]> nhitcut[det1][i1]) fGainFact[det1][i1][j1][k1]=-99.0;
373                    if(nhitcell[det1][i1][j1][k1]< nhitcut[det1][i1])
374                      {
375                       isoMean[det1][i1][j1][k1]=fHadcIso[det1][i1][j1][k1]->
376                         GetMean();
377                       if(isoMean[det1][i1][j1][k1]>0.0 && histMean[det1][i1]>0.0)
378                         {
379                           fGainFact[det1][i1][j1][k1]=
380                             isoMean[det1][i1][j1][k1]/(histMean[det1][i1]*
381                                                      smNormFactor[det1][i1]);
382
383                         }
384                     }   
385                    Float_t gain = fGainFact[det1][i1][j1][k1];
386                    fCalibGain->SetGainFact(det1,i1,j1,k1,gain);
387                 }
388                           
389                         
390             }
391         }
392     }
393 }//CalculateIsoCell()
394
395 // ------------------------------------------------------------------------ //
396
397 Bool_t AliPMDCalibrator::Store()
398 {
399   AliCDBManager *man = AliCDBManager::Instance();
400   //man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
401   if(!man->IsDefaultStorageSet()) return kFALSE;
402   AliCDBId id("PMD/Calib/Gain",0,999999999);
403   AliCDBMetaData md;
404   md.SetBeamPeriod(0);
405   md.SetComment("Test");
406   
407   printf("\n\n\n fCalibData\n");
408   //fCalibData->Print(0);
409   //printf("\n\n\n fCalibData\n");
410   
411   Bool_t result = man->Put(fCalibGain,id,&md);
412
413   return result;
414 }
415