added method to read pedestal
[u/mrichter/AliRoot.git] / PMD / AliPMDCalibGain.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 "TTree.h"
27
28 // --- Standard library ---
29
30 // --- AliRoot header files ---
31 #include "AliDAQ.h"
32 #include "AliLog.h"
33 #include "AliRawReader.h"
34 #include "AliPMDRawStream.h"
35 #include "AliPMDddldata.h"
36 #include "AliPMDCalibGain.h"
37
38 ClassImp(AliPMDCalibGain)
39
40 AliPMDCalibGain::AliPMDCalibGain():
41   TObject(),
42   fpw(NULL)
43 {
44   // Standard Constructor
45     for(Int_t idet = 0; idet < kDet; idet++)
46     {
47      fDetCount[kDet] =0.;
48         for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
49         {
50             fSMIso[idet][ismn]   = 0.;
51             fSMCount[idet][ismn] = 0.;
52             for(Int_t jrow = 0; jrow < kMaxRow; jrow++)
53             {
54                 for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
55                 {
56                     fCellIso[idet][ismn][jrow][kcol]   = 0.;
57                     fCellCount[idet][ismn][jrow][kcol] = 0.;
58                     fPedMeanRMS[idet][ismn][jrow][kcol] = 0.;
59
60                 }
61             }
62         }
63     }
64
65
66
67 }
68 // ------------------------------------------------------------------------ //
69 AliPMDCalibGain::AliPMDCalibGain(const AliPMDCalibGain &pmdcalibgain):
70   TObject(pmdcalibgain),
71   fpw(NULL)
72 {
73     for(Int_t idet = 0; idet < kDet; idet++)
74     {
75      fDetCount[idet] = pmdcalibgain.fDetCount[idet];
76      fDetIso[idet] = pmdcalibgain.fDetIso[idet];
77         for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
78         {
79             fSMIso[idet][ismn] = pmdcalibgain.fSMIso[idet][ismn] ;
80             fSMCount[idet][ismn] = pmdcalibgain.fSMCount[idet][ismn] ;
81             for(Int_t jrow = 0; jrow < kMaxRow; jrow++)
82             {
83                 for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
84                 {
85                   fCellIso[idet][ismn][jrow][kcol]    =
86                     pmdcalibgain.fCellIso[idet][ismn][jrow][kcol];
87                   fCellCount[idet][ismn][jrow][kcol]  =
88                     pmdcalibgain.fCellCount[idet][ismn][jrow][kcol];
89                   fPedMeanRMS[idet][ismn][jrow][kcol] =
90                     pmdcalibgain.fPedMeanRMS[idet][ismn][jrow][kcol];
91
92                 }
93             }
94         }
95     }
96     
97 }
98 // ------------------------------------------------------------------------ //
99 AliPMDCalibGain &AliPMDCalibGain::operator=(const AliPMDCalibGain &pmdcalibgain)
100 {
101     if(this != &pmdcalibgain)
102     {
103       this->fpw = pmdcalibgain.fpw;
104         for(Int_t idet = 0; idet < kDet; idet++)
105         {
106             for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
107             {
108                 fSMIso[idet][ismn] = pmdcalibgain.fSMIso[idet][ismn];
109                 fSMCount[idet][ismn] = pmdcalibgain.fSMCount[idet][ismn];
110                 for(Int_t jrow = 0; jrow < kMaxRow;jrow++)
111                 {
112                     for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
113                       {
114                         fCellIso[idet][ismn][jrow][kcol]  =
115                           pmdcalibgain.fCellIso[idet][ismn][jrow][kcol];
116                         fCellCount[idet][ismn][jrow][kcol]  =
117                           pmdcalibgain.fCellCount[idet][ismn][jrow][kcol];
118                         fPedMeanRMS[idet][ismn][jrow][kcol] = 
119                           pmdcalibgain.fPedMeanRMS[idet][ismn][jrow][kcol];
120
121                     }
122                 }
123             }
124         }
125     }
126     return *this;
127 }
128 // ------------------------------------------------------------------------ //
129 AliPMDCalibGain::~AliPMDCalibGain()
130 {
131     // dtor
132
133 }
134
135 // ------------------------------------------------------------------------ //
136
137 Int_t AliPMDCalibGain::ExtractPedestal()
138 {
139   // Pedestal extraction from the PMD_PED.root file
140   // To be called once at the beginning
141
142   Int_t   det, sm, row, col;
143   Float_t mean, rms;
144
145   TFile *pedfile = new TFile("PMD_PED1.root");
146
147   if(!pedfile)
148     {
149       printf("ERROR --- NO PEDESTAL (PMD_PED1.root) FILE IS FOUND --- STOP GAIN DA\n");
150       return -3;
151     }
152
153
154   TTree *ped =(TTree*)pedfile->Get("ped");
155
156   ped->SetBranchAddress("det",&det);
157   ped->SetBranchAddress("sm",&sm);
158   ped->SetBranchAddress("row",&row);
159   ped->SetBranchAddress("col",&col);
160   ped->SetBranchAddress("mean",&mean);
161   ped->SetBranchAddress("rms",&rms);
162
163   Int_t nentries = (Int_t)ped->GetEntries();
164
165   for (Int_t ient = 0; ient < nentries; ient++)
166     {
167       ped->GetEntry(ient);
168       fPedMeanRMS[det][sm][row][col] = mean + 3.*rms;
169       //printf("Mean= %f, RMS= %f, PedMeanRMS=%f\n",mean,rms,fPedMeanRMS[det][sm][row][col]);
170
171     }
172
173   pedfile->Close();
174   delete pedfile;
175   pedfile = 0x0;
176
177   return 1;
178 }
179 // ------------------------------------------------------------------------ //
180
181 void AliPMDCalibGain::ReadIntermediateFile()
182 {
183   // Read the variables from the file
184   
185   fpw = fopen("interfile.dat","r");
186
187   Float_t detcount, detiso;
188   Float_t smcount, smiso;
189   Float_t cellcount, celliso;
190
191   for (Int_t idet = 0; idet < kDet; idet++)
192     {
193       fscanf(fpw,"%d %f %f",&idet,&detcount,&detiso);
194       fDetCount[idet] = detcount;
195       fDetIso[idet]   = detiso;
196     }
197
198   for (Int_t idet = 0; idet < kDet; idet++)
199     {
200       for (Int_t ism = 0; ism < kMaxSMN; ism++)
201         {
202           fscanf(fpw,"%d %d %f %f",&idet,&ism,&smcount,&smiso);
203
204           fSMCount[idet][ism] = smcount;
205           fSMIso[idet][ism]   = smiso;
206         }
207     }
208
209   for (Int_t idet = 0; idet < kDet; idet++)
210     {
211       for (Int_t ism = 0; ism < kMaxSMN; ism++)
212         {
213           for (Int_t irow = 0; irow < kMaxRow; irow++)
214             {
215               for (Int_t icol = 0; icol < kMaxCol; icol++)
216                 {
217                   fscanf(fpw,"%d %d %d %d %f %f",&idet,&ism,&irow,&icol,
218                           &cellcount,&celliso);
219
220                   fCellCount[idet][ism][irow][icol] = cellcount;
221                   fCellIso[idet][ism][irow][icol]   = celliso;
222                 }
223             }
224         }
225     }
226
227   fclose(fpw);
228
229 }
230 // ------------------------------------------------------------------------ //
231 void AliPMDCalibGain::WriteIntermediateFile()
232 {
233
234  //Following variables to be written
235
236   /*
237     fDetIso[idet] ;
238     fSMIso[idet][ismn]; 
239     fCellIso[idet][ismn][irow][icol]; 
240     
241     fDetCount[idet];
242     fSMCount[idet][ismn];
243     fCellCount[idet][ismn][irow][icol];
244   */                              
245
246   fpw = fopen("interfile.dat","w+");
247
248   for (Int_t idet = 0; idet < kDet; idet++)
249     {
250       fprintf(fpw,"%d %f %f\n",idet,fDetCount[idet],fDetIso[idet]);
251     }
252
253   for (Int_t idet = 0; idet < kDet; idet++)
254     {
255       for (Int_t ism = 0; ism < kMaxSMN; ism++)
256         {
257           fprintf(fpw,"%d %d %f %f\n",idet,ism, fSMCount[idet][ism],fSMIso[idet][ism]);
258         }
259     }
260
261   for (Int_t idet = 0; idet < kDet; idet++)
262     {
263       for (Int_t ism = 0; ism < kMaxSMN; ism++)
264         {
265           for (Int_t irow = 0; irow < kMaxRow; irow++)
266             {
267               for (Int_t icol = 0; icol < kMaxCol; icol++)
268                 {
269                   fprintf(fpw,"%d %d %d %d %f %f\n",idet,ism,irow,icol,
270                           fCellCount[idet][ism][irow][icol],
271                           fCellIso[idet][ism][irow][icol]);
272                 }
273             }
274         }
275     }
276
277   fclose(fpw);
278
279 }
280
281 // ------------------------------------------------------------------------ //
282
283 Bool_t AliPMDCalibGain::ProcessEvent(AliRawReader *rawReader, TObjArray *pmdddlcont)
284 {
285   // Calculates the ADC of isolated cell
286
287   const Int_t kDDL           = AliDAQ::NumberOfDdls("PMD");
288   const Int_t kCellNeighbour = 6;
289   Int_t neibx[6] = {1,0,-1,-1,0,1};
290   Int_t neiby[6] = {0,1,1,0,-1,-1};
291   
292   Int_t id1,jd1;  //neighbour row/col
293   Int_t isocount; //number of neighbours with 0 signal
294
295   Float_t d1[kDet][kMaxSMN][kMaxRow][kMaxCol];
296
297   for(Int_t idet = 0; idet < kDet; idet++)
298     {
299       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
300         {
301           for(Int_t irow = 0; irow < kMaxRow; irow++)
302             {
303               for(Int_t icol = 0; icol < kMaxCol; icol++)
304                 {
305                   d1[idet][ismn][irow][icol] = 0.;
306                 }
307             }
308         }
309     }
310
311   AliPMDRawStream rawStream(rawReader);
312
313   Int_t iddl = -1;
314
315   Int_t numberofDDLs = 0;
316
317     while ((iddl = rawStream.DdlData(pmdddlcont)) >=0) {
318       numberofDDLs++;
319
320       Int_t ientries = pmdddlcont->GetEntries();
321
322       for (Int_t ient = 0; ient < ientries; ient++)
323       {
324           AliPMDddldata *pmdddl = (AliPMDddldata*)pmdddlcont->UncheckedAt(ient);
325           
326           Int_t idet = pmdddl->GetDetector();
327           Int_t ismn = pmdddl->GetSMN();
328           //Int_t mcm = pmdddl->GetMCM();
329           //Int_t ichno = pmdddl->GetChannel();
330           Int_t irow = pmdddl->GetRow();
331           Int_t icol = pmdddl->GetColumn();
332           Int_t isig = pmdddl->GetSignal();
333
334
335           // Pedestal subtraction
336
337           if (isig>0)
338             {
339               d1[idet][ismn][irow][icol] =
340                 (Float_t) isig - fPedMeanRMS[idet][ismn][irow][icol];
341 //printf("Signal_ped_subtracted=%f, pedestal=%f\n",d1[idet][ismn][irow][icol]),fPedMeanRMS[idet][ismn][irow][icol];
342             }
343       }
344       pmdddlcont->Delete();
345   }
346   
347   for(Int_t idet=0; idet < kDet; idet++)
348   {
349       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
350       {
351           for(Int_t irow = 0; irow < kMaxRow; irow++)
352           {
353               for(Int_t icol = 0; icol < kMaxCol; icol++)
354               {
355                   if(d1[idet][ismn][irow][icol] > 0)
356                   {
357                       isocount = 0;
358                       for(Int_t ii = 0; ii < kCellNeighbour; ii++)
359                       {
360                           id1 = irow + neibx[ii];
361                           jd1 = icol + neiby[ii];
362                           if(d1[idet][ismn][id1][jd1] == 0)
363                           {
364                               isocount++;
365                               if(isocount == kCellNeighbour)
366                               {
367                                   fDetIso[idet] += d1[idet][ismn][irow][icol];
368                                   fSMIso[idet][ismn] += d1[idet][ismn][irow][icol];
369                                   fCellIso[idet][ismn][irow][icol] += d1[idet][ismn][irow][icol];
370                                   fDetCount[idet]++;
371                                   fSMCount[idet][ismn]++;
372                                   fCellCount[idet][ismn][irow][icol]++;
373                                   
374                               }
375                           }
376                       }  // neigh cell cond.
377                   }     // d>0 cond.
378               }
379           }
380       }
381   }
382
383   if (numberofDDLs < kDDL)
384       return kFALSE;
385   return kTRUE;
386
387 }
388 // ------------------------------------------------------------------------ //
389 void AliPMDCalibGain::Analyse(TTree *gaintree)
390 {
391     // Calculates the mean
392     Int_t   det, sm, row, col;
393     Float_t gain;
394     Float_t modmean  = 0.;
395     Float_t cellmean = 0.;
396     Float_t detmean =0.;
397
398     gaintree->Branch("det",&det,"det/I");
399     gaintree->Branch("sm",&sm,"sm/I");
400     gaintree->Branch("row",&row,"row/I");
401     gaintree->Branch("col",&col,"col/I");
402     gaintree->Branch("gain",&gain,"gain/F");
403
404     for(Int_t idet = 0; idet < kDet; idet++)
405       {
406         if (fDetCount[idet]>0 )
407           detmean=fDetIso[idet]/fDetCount[idet];
408         for(Int_t ism = 0; ism < kMaxSMN; ism++)
409           {
410             if (fSMCount[idet][ism] > 0)
411               modmean = fSMIso[idet][ism]/fSMCount[idet][ism];
412             for(Int_t irow = 0; irow < kMaxRow; irow++)
413               {
414                 for(Int_t icol = 0; icol < kMaxCol; icol++)
415                   {
416                     
417                     if (fCellCount[idet][ism][irow][icol] > 0.)
418                       {
419                         cellmean = fCellIso[idet][ism][irow][icol]/fCellCount[idet][ism][irow][icol];
420                       }
421                     det      = idet;
422                     sm       = ism;
423                     row      = irow;
424                     col      = icol;
425                     if (cellmean > 0.0 && fCellCount[idet][ism][irow][icol]>0.)
426                       {
427                         gain = cellmean/detmean;
428                       }
429                     else
430                       {
431                         gain = -1.;
432                       }
433                     //if(fCellCount[idet][ism][irow][icol]>0.) printf("CellCount =%f, gain= %f\n",fCellCount[idet][ism][irow][icol],gain);
434                     gaintree->Fill();
435                   }
436               }
437           }
438       }
439     
440 }
441 // ------------------------------------------------------------------------ //