warning fixed
[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   if(fHdetIso) delete fHdetIso ;
127   if(fHsmIso)  delete fHsmIso ;
128   if(fHadcIso) delete fHadcIso ;
129   delete fCalibGain;
130   delete fCalibPed;
131 }
132 // ------------------------------------------------------------------------ //
133
134 void AliPMDCalibrator::Exec()
135 {
136   // reads parameters and does the calibration
137   CalculateIsoCell() ;
138 }
139 // ------------------------------------------------------------------------ //
140
141 void AliPMDCalibrator::Init()
142 {
143   // intializes everything
144   char hname2[kDet];
145   char htitle2[2];
146   char hname[kMaxSMN];
147   char hname24[kMaxSMN];
148   char hnameiso[120];
149   char htitle1[120];
150   
151   for(Int_t d = 0; d < kDet; d++) {
152     sprintf(hname2,"Isolated cell adc for Det Plane %d",d);
153     fHdetIso[d]= new TH1F(hname2,htitle2,100,0,1000);
154     for(Int_t i1 = 0; i1 < kMaxSMN; i1++) {
155       sprintf(hname,"det_%d_iso_sm_%2d",d,i1);
156       sprintf(hname24,"det_%d_iso_sm_%2d",d,i1);
157       fHsmIso[d][i1]= new TH1F(hname,hname24,100,0,1000);
158       for(Int_t j1 = 0; j1 < kMaxRow; j1++) {
159         for(Int_t k1 = 0; k1 < kMaxCol; k1++) {
160           sprintf(hnameiso,"Isolated Cell ADC for det_%d_cell_sm%d_row%d_col%d"
161                   ,d,i1,j1,k1);
162           sprintf(htitle1,"Isolated Cell ADC for det_%d_cell_sm%d_row%d_col%d"
163                   ,d,i1,j1,k1);
164           
165           TObject *old=gDirectory->GetList()->FindObject(hnameiso);
166           if (old) gDirectory->GetList()->Remove(old);
167           fHadcIso[d][i1][j1][k1] = new TH1F(hnameiso,htitle1,100,0.,1000.);
168         }
169       }
170     }
171   }
172   
173 }
174
175 // ------------------------------------------------------------------------ //
176
177 void AliPMDCalibrator::CalculateIsoCell()
178 {
179   // Calculates the ADC of isolated cell
180
181   TObjArray pmdddlcont;
182
183   const Int_t kCellNeighbour = 6;
184
185   Int_t neibx[6] = {1,0,-1,-1,0,1};
186   Int_t neiby[6] = {0,1,1,0,-1,-1};
187
188   Int_t id1,jd1;            //neighbour row/col
189   Int_t countisocell = 0 ;  //number of isilated cell
190   Int_t isocount;           //number of neighbours with 0 signal
191   Int_t d1[kDet][kMaxSMN][kMaxRow][kMaxCol];
192   Int_t maxhit;
193   Int_t nhit[kDet][kMaxSMN];
194   Int_t nhitcell[kDet][kMaxSMN][kMaxRow][kMaxCol];
195   
196   for(Int_t idet = 0; idet < kDet; idet++)
197     {
198       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
199         {
200           nhit[idet][ismn] = 0;
201           for(Int_t irow = 0; irow < kMaxRow; irow++)
202             {
203               for(Int_t icol = 0; icol < kMaxCol; icol++)
204                 {
205                   d1[idet][ismn][irow][icol] = 0;
206                   nhitcell[idet][ismn][irow][icol] = 0;
207                 }
208             }
209         }
210     }
211
212   Float_t tempnhit1[kDet][kMaxSMN];
213   Float_t tempnhit2[kDet][kMaxSMN];
214   Float_t meannhit[kDet][kMaxSMN];
215   Float_t meanSqnhit[kDet][kMaxSMN];
216   Float_t sigmanhit[kDet][kMaxSMN];
217   Float_t count[kDet][kMaxSMN];
218   Float_t nhitcut[kDet][kMaxSMN];
219
220   for (Int_t idet = 0; idet < kDet; idet++)
221   {
222       for (Int_t ismn = 0; ismn < kMaxSMN; ismn++)
223       {
224           tempnhit1[idet][ismn]  = 0.;
225           tempnhit2[idet][ismn]  = 0.;
226           meannhit[idet][ismn]   = 0.;
227           meanSqnhit[idet][ismn] = 0.;
228           sigmanhit[idet][ismn]  = 0.;
229           count[idet][ismn]      = 0.;
230           nhitcut[idet][ismn]    = 0.;
231       }
232   }
233
234
235   //accessing raw data
236   AliRawReaderFile reader(".");
237   AliPMDRawStream stream(&reader);
238   while(reader.NextEvent())
239     { 
240       // New PMD Reader is plugged in
241         Int_t iddl = -1;
242         while ((iddl = stream.DdlData(&pmdddlcont)) >=0) {
243
244             //reader.Select("PMD", iddl, iddl);
245             //stream.DdlData(iddl,&pmdddlcont);
246           Int_t ientries = pmdddlcont.GetEntries();
247           for (Int_t ient = 0; ient < ientries; ient++)
248             {
249               AliPMDddldata *pmdddl = 
250                 (AliPMDddldata*)pmdddlcont.UncheckedAt(ient);
251               Int_t idet = pmdddl->GetDetector();
252               Int_t ismn = pmdddl->GetSMN();
253               Int_t irow = pmdddl->GetRow();
254               Int_t icol = pmdddl->GetColumn();
255               Float_t isig1 = pmdddl->GetSignal();
256               // Pedestal Subtraction
257               Int_t   pedmeanrms = 
258                 fCalibPed->GetPedMeanRms(idet,ismn,irow,icol);
259               Int_t   pedrms1    = (Int_t) pedmeanrms%100;
260               Float_t pedrms     = (Float_t)pedrms1/10.;
261               Float_t pedmean    = (Float_t) (pedmeanrms - pedrms1)/1000.0;
262               Float_t isig = isig1 - (pedmean + 3.0*pedrms);
263               if (isig>0)
264                 {
265                   d1[idet][ismn][irow][icol] = (Int_t)isig;
266                   nhitcell[idet][ismn][irow][icol] += 1; 
267                 }
268             }//ient loop
269           pmdddlcont.Clear();
270         }//iddl loop
271       maxhit = 0;
272       for(Int_t idet=0; idet < kDet; idet++)
273         {
274           for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
275             {
276               for(Int_t irow = 0; irow < kMaxRow; irow++)
277                 {
278                   for(Int_t icol = 0; icol < kMaxCol; icol++)
279                     {
280                       if(d1[idet][ismn][irow][icol] > 0)
281                         {
282                           isocount = 0;
283                           for(Int_t ii = 0; ii < kCellNeighbour; ii++)
284                             {
285                               id1 = irow + neibx[ii];
286                               jd1 = icol + neiby[ii];
287                               if(d1[idet][ismn][id1][jd1] == 0)
288                                 {
289                                   isocount++;
290                                   if(isocount == kCellNeighbour)
291                                     {
292                                       countisocell++;
293                                       maxhit++;
294                                       fHdetIso[idet]->
295                                         Fill(d1[idet][ismn][irow][icol]);
296                                       fHsmIso[idet][ismn]->
297                                         Fill(d1[idet][ismn][irow][icol]);
298                                       fHadcIso[idet][ismn][irow][icol]->
299                                         Fill(d1[idet][ismn][irow][icol]);
300                                     }
301                                 }
302                             }  // neigh cell cond.
303                         }     // d>0 cond.
304                     }
305                 }
306             }
307         } //det
308     }//event loop
309
310   // Mean and Sigma Calculations
311   for(Int_t idet=0; idet < kDet; idet++){
312       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++){
313           for(Int_t irow = 0; irow < kMaxRow; irow++){
314               for(Int_t icol = 0; icol < kMaxCol; icol++){
315                   if(nhitcell[idet][ismn][irow][icol]>0){
316                       count[idet][ismn] += 1;
317                       tempnhit1[idet][ismn] += nhitcell[idet][ismn][irow][icol];
318                       tempnhit2[idet][ismn] += nhitcell[idet][ismn][irow][icol]
319                           *nhitcell[idet][ismn][irow][icol];
320                   }
321               }
322           }
323       }
324   }//det loop
325
326   //cout<<"nhit cell = "<<idet<<"  "<<ismn<<"  "
327   //  <<irow<<"  "<<icol<<"  "<<nhitcell[idet][ismn][irow][icol]<<endl;
328   //count[idet][ismn] += 1;
329
330   for(Int_t i=0; i < kDet; i++)
331     {
332       for(Int_t j=0; j < kMaxSMN; j++)
333         {
334           if(count[i][j] > 0.)
335             {
336                 meannhit[i][j]   = tempnhit1[i][j]/count[i][j];
337                 meanSqnhit[i][j] = tempnhit2[i][j]/count[i][j]; 
338                 sigmanhit[i][j]  = sqrt(meanSqnhit[i][j]-
339                                       (meannhit[i][j]*meannhit[i][j]));
340                 nhitcut[i][j]    = 3*sigmanhit[i][j] + meannhit[i][j];
341             }
342         }
343     }
344
345
346   Double_t histdetMean[kDet];
347   Double_t histMean[kDet][kMaxSMN];
348   Double_t isoMean[kDet][kMaxSMN][kMaxRow][kMaxCol];
349   Double_t smNormFactor[kDet][kMaxSMN];
350   
351   for(Int_t det1 = 0; det1 < kDet; det1++)
352     {
353       histdetMean[det1]= fHdetIso[det1]->GetMean();
354       for(Int_t i1 = 0; i1 < kMaxSMN; i1++)
355         {
356           histMean[det1][i1]= fHsmIso[det1][i1]->GetMean();
357           if(histMean[det1][i1]>0.0 && histdetMean[det1]>0.0)
358             {
359               smNormFactor[det1][i1]= histdetMean[det1]/histMean[det1][i1];
360             }
361           for(Int_t j1 = 0; j1 < kMaxRow; j1++)
362             {
363               for(Int_t k1 = 0; k1 < kMaxCol; k1++)
364                 {
365                  if(nhitcell[det1][i1][j1][k1]> nhitcut[det1][i1]) fGainFact[det1][i1][j1][k1]=-99.0;
366                    if(nhitcell[det1][i1][j1][k1]< nhitcut[det1][i1])
367                      {
368                       isoMean[det1][i1][j1][k1]=fHadcIso[det1][i1][j1][k1]->
369                         GetMean();
370                       if(isoMean[det1][i1][j1][k1]>0.0 && histMean[det1][i1]>0.0)
371                         {
372                           fGainFact[det1][i1][j1][k1]=
373                             isoMean[det1][i1][j1][k1]/(histMean[det1][i1]*
374                                                      smNormFactor[det1][i1]);
375
376                         }
377                     }   
378                    Float_t gain = fGainFact[det1][i1][j1][k1];
379                    fCalibGain->SetGainFact(det1,i1,j1,k1,gain);
380                 }
381                           
382                         
383             }
384         }
385     }
386 }//CalculateIsoCell()
387
388 // ------------------------------------------------------------------------ //
389
390 Bool_t AliPMDCalibrator::Store()
391 {
392   AliCDBManager *man = AliCDBManager::Instance();
393   //man->SetDefaultStorage("local://$ALICE_ROOT");
394   if(!man->IsDefaultStorageSet()) return kFALSE;
395   AliCDBId id("PMD/Calib/Gain",0,999999999);
396   AliCDBMetaData md;
397   md.SetBeamPeriod(0);
398   md.SetComment("Test");
399   
400   printf("\n\n\n fCalibData\n");
401   //fCalibData->Print(0);
402   //printf("\n\n\n fCalibData\n");
403   
404   Bool_t result = man->Put(fCalibGain,id,&md);
405
406   return result;
407 }
408