Bug fixed to take care of thearray length going beyond the limit
[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 #include "TMath.h"
28
29 // --- Standard library ---
30
31 // --- AliRoot header files ---
32 #include "AliDAQ.h"
33 #include "AliLog.h"
34 #include "AliRawReader.h"
35 #include "AliPMDRawStream.h"
36 #include "AliPMDddldata.h"
37 #include "AliPMDCalibGain.h"
38
39 ClassImp(AliPMDCalibGain)
40
41 AliPMDCalibGain::AliPMDCalibGain():
42   TObject(),
43   fpw(NULL)
44 {
45   // Standard Constructor
46   for(Int_t idet = 0; idet < kDet; idet++)
47     {
48       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
49         {
50           fSMIso[idet][ismn]   = 0.;
51           fSMCount[idet][ismn] = 0.;
52           fCountSm[idet][ismn]=0.;
53           fTempnhit[idet][ismn]=0.;
54           fTempnhitSq[idet][ismn]=0.;
55           for(Int_t jrow = 0; jrow < kMaxRow; jrow++)
56             {
57               for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
58                 {
59                   fCellIso[idet][ismn][jrow][kcol]    = 0.;
60                   fCellCount[idet][ismn][jrow][kcol]  = 0.;
61                   fNhitCell[idet][ismn][jrow][kcol]   = 0.;
62                   fPedMeanRMS[idet][ismn][jrow][kcol] = 0.;
63                   fHotFlag[idet][ismn][jrow][kcol]    = 0.;
64                   
65                 }
66             }
67         }
68     }
69   
70 }
71 // ------------------------------------------------------------------------ //
72 AliPMDCalibGain::AliPMDCalibGain(const AliPMDCalibGain &pmdcalibgain):
73   TObject(pmdcalibgain),
74   fpw(NULL)
75 {
76   for(Int_t idet = 0; idet < kDet; idet++)
77     {
78       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
79         {
80           fSMIso[idet][ismn] = pmdcalibgain.fSMIso[idet][ismn] ;
81           fSMCount[idet][ismn] = pmdcalibgain.fSMCount[idet][ismn] ;
82           fCountSm[idet][ismn] = pmdcalibgain.fCountSm[idet][ismn];
83           fTempnhit[idet][ismn] = pmdcalibgain.fTempnhit[idet][ismn];
84           fTempnhitSq[idet][ismn] = pmdcalibgain.fTempnhitSq[idet][ismn];
85           for(Int_t jrow = 0; jrow < kMaxRow; jrow++)
86             {
87               for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
88                 {
89                   fCellIso[idet][ismn][jrow][kcol] = pmdcalibgain.fCellIso[idet][ismn][jrow][kcol];
90                   fCellCount[idet][ismn][jrow][kcol] = pmdcalibgain.fCellCount[idet][ismn][jrow][kcol];
91                   fNhitCell[idet][ismn][jrow][kcol] = pmdcalibgain.fNhitCell[idet][ismn][jrow][kcol];
92                   fPedMeanRMS[idet][ismn][jrow][kcol] = pmdcalibgain.fPedMeanRMS[idet][ismn][jrow][kcol];
93                   fHotFlag[idet][ismn][jrow][kcol] = pmdcalibgain.fHotFlag[idet][ismn][jrow][kcol];
94                   
95                 }
96             }
97         }
98     }
99   
100 }
101 // ------------------------------------------------------------------------ //
102 AliPMDCalibGain &AliPMDCalibGain::operator=(const AliPMDCalibGain &pmdcalibgain)
103 {
104   if(this != &pmdcalibgain)
105     {
106       this->fpw = pmdcalibgain.fpw;
107       for(Int_t idet = 0; idet < kDet; idet++)
108         {
109           for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
110             {
111               fSMIso[idet][ismn]      = pmdcalibgain.fSMIso[idet][ismn];
112               fSMCount[idet][ismn]    = pmdcalibgain.fSMCount[idet][ismn];
113               fCountSm[idet][ismn]    = pmdcalibgain.fCountSm[idet][ismn];
114               fTempnhit[idet][ismn]   = pmdcalibgain.fTempnhit[idet][ismn];
115               fTempnhitSq[idet][ismn] = pmdcalibgain.fTempnhitSq[idet][ismn];
116               for(Int_t jrow = 0; jrow < kMaxRow;jrow++)
117                 {
118                   for(Int_t kcol = 0; kcol < kMaxCol; kcol++)
119                     {
120                       fCellIso[idet][ismn][jrow][kcol] = pmdcalibgain.fCellIso[idet][ismn][jrow][kcol];
121                       fCellCount[idet][ismn][jrow][kcol] = pmdcalibgain.fCellCount[idet][ismn][jrow][kcol];
122                       fNhitCell[idet][ismn][jrow][kcol] = pmdcalibgain.fNhitCell[idet][ismn][jrow][kcol]; //za 
123                       fPedMeanRMS[idet][ismn][jrow][kcol] = pmdcalibgain.fPedMeanRMS[idet][ismn][jrow][kcol];
124                       fHotFlag[idet][ismn][jrow][kcol] = pmdcalibgain.fHotFlag[idet][ismn][jrow][kcol];
125                       
126                     }
127                 }
128             }
129         }
130     }
131   return *this;
132 }
133 // ------------------------------------------------------------------------ //
134 AliPMDCalibGain::~AliPMDCalibGain()
135 {
136   // dtor
137 }
138
139 // ------------------------------------------------------------------------ //
140
141 Int_t AliPMDCalibGain::ExtractPedestal(const Char_t *rootFile)
142 {
143   // Pedestal extraction from the PMD_PED.root file
144   // To be called once at the beginning
145   
146   Int_t   det, sm, row, col;
147   Float_t mean, rms;
148   
149   TFile *pedfile = new TFile(rootFile);
150   
151   if(!pedfile)
152     {
153       printf("ERROR --- NO PEDESTAL (PMD_PED1.root) FILE IS FOUND --- STOP GAIN DA\n");
154       return -3;
155     }
156
157
158   TTree *ped =(TTree*)pedfile->Get("ped");
159
160   ped->SetBranchAddress("det",&det);
161   ped->SetBranchAddress("sm",&sm);
162   ped->SetBranchAddress("row",&row);
163   ped->SetBranchAddress("col",&col);
164   ped->SetBranchAddress("mean",&mean);
165   ped->SetBranchAddress("rms",&rms);
166
167   Int_t nentries = (Int_t)ped->GetEntries();
168
169   for (Int_t ient = 0; ient < nentries; ient++)
170     {
171       ped->GetEntry(ient);
172       fPedMeanRMS[det][sm][row][col] = mean + 3.*rms;
173       //printf("Mean= %f, RMS= %f, PedMeanRMS=%f\n",mean,rms,fPedMeanRMS[det][sm][row][col]);
174
175     }
176
177   pedfile->Close();
178   delete pedfile;
179   pedfile = 0x0;
180
181   return 1;
182 }
183 //------------------------------------------------------------------------------------------------
184
185 Int_t AliPMDCalibGain::ExtractHotChannel(const Char_t *rootFile)
186 {
187   // HotChannel extraction from the PMD_HOT.root file
188   // To be called once at the beginning
189
190   Int_t   det, sm, row, col;
191   Float_t flag;
192
193   TFile *hotmapfile = new TFile(rootFile);
194
195   if(!hotmapfile)
196     {
197       printf(" NO HOTCHANNEL MAP FOUND (PMD_HOT.root) FILE IS FOUND \n");
198       fHotFlag[kDet][kMaxSMN][kMaxRow][kMaxCol] = 0.;
199     }
200
201
202   TTree *hot =(TTree*)hotmapfile->Get("hot");
203   
204   hot->SetBranchAddress("det",&det);
205   hot->SetBranchAddress("sm",&sm);
206   hot->SetBranchAddress("row",&row);
207   hot->SetBranchAddress("col",&col);
208   hot->SetBranchAddress("flag",&flag);
209   
210   Int_t nentries = (Int_t)hot->GetEntries();
211   
212   for (Int_t ient = 0; ient < nentries; ient++)
213     {
214       hot->GetEntry(ient);
215       fHotFlag[det][sm][row][col] = flag;
216
217       //printf(" HotFlag=%f\n",fHotFlag[det][sm][row][col]);
218     }
219
220   hotmapfile->Close();
221   delete hotmapfile;
222   hotmapfile = 0x0;
223   
224   return 1;
225 }
226
227
228 // ------------------------------------------------------------------------ //
229
230 void AliPMDCalibGain::ReadTempFile(const Char_t *tempFile)
231 {
232   // Read the variables from the file
233   
234   fpw = fopen(tempFile,"r");
235   
236   Float_t smcount, smiso;
237   Float_t cellcount, celliso;
238
239
240   for (Int_t idet = 0; idet < kDet; idet++)
241     {
242       for (Int_t ism = 0; ism < kMaxSMN; ism++)
243         {
244           fscanf(fpw,"%d %d %f %f",&idet,&ism,&smcount,&smiso);
245           
246           fSMCount[idet][ism] = smcount;
247           fSMIso[idet][ism]   = smiso;
248         }
249     }
250
251   for (Int_t idet = 0; idet < kDet; idet++)
252     {
253       for (Int_t ism = 0; ism < kMaxSMN; ism++)
254         {
255           for (Int_t irow = 0; irow < kMaxRow; irow++)
256             {
257               for (Int_t icol = 0; icol < kMaxCol; icol++)
258                 {
259                   fscanf(fpw,"%d %d %d %d %f %f",&idet,&ism,&irow,&icol,
260                          &cellcount,&celliso);
261                   
262                   fCellCount[idet][ism][irow][icol] = cellcount;
263                   fCellIso[idet][ism][irow][icol]   = celliso;
264                 }
265             }
266         }
267     }
268   
269   fclose(fpw);
270
271 }
272 // ------------------------------------------------------------------------ //
273 void AliPMDCalibGain::WriteTempFile(const Char_t *tempFile)
274 {
275   // Write the Temporary file if the required statics is not achieved
276
277
278   /*
279     Following variables to be written to a file
280     fDetIso[idet] ;
281     fSMIso[idet][ismn]; 
282     fCellIso[idet][ismn][irow][icol]; 
283     
284     fDetCount[idet];
285     fSMCount[idet][ismn];
286     fCellCount[idet][ismn][irow][icol];
287   */                              
288
289
290   fpw = fopen(tempFile,"w+");
291   
292   for (Int_t idet = 0; idet < kDet; idet++)
293     {
294      // fprintf(fpw,"%d %f %f\n",idet,fDetCount[idet],fDetIso[idet]);
295     }
296   
297   for (Int_t idet = 0; idet < kDet; idet++)
298     {
299       for (Int_t ism = 0; ism < kMaxSMN; ism++)
300         {
301           fprintf(fpw,"%d %d %f %f\n",idet,ism, fSMCount[idet][ism],fSMIso[idet][ism]);
302         }
303     }
304   
305   for (Int_t idet = 0; idet < kDet; idet++)
306     {
307       for (Int_t ism = 0; ism < kMaxSMN; ism++)
308         {
309           for (Int_t irow = 0; irow < kMaxRow; irow++)
310             {
311               for (Int_t icol = 0; icol < kMaxCol; icol++)
312                 {
313                   fprintf(fpw,"%d %d %d %d %f %f\n",idet,ism,irow,icol,
314                           fCellCount[idet][ism][irow][icol],
315                           fCellIso[idet][ism][irow][icol]);
316                 }
317             }
318         }
319     }
320   
321   fclose(fpw);
322   
323 }
324
325 // ------------------------------------------------------------------------ //
326
327 Bool_t AliPMDCalibGain::ProcessEvent(AliRawReader *rawReader, TObjArray *pmdddlcont)
328 {
329   // Calculates the ADC of isolated cell
330   
331   const Int_t kDDL           = AliDAQ::NumberOfDdls("PMD");
332   const Int_t kCellNeighbour = 6;
333   Int_t neibx[6] = {1,0,-1,-1,0,1};
334   Int_t neiby[6] = {0,1,1,0,-1,-1};
335   
336   Int_t id1,jd1;  //neighbour row/col
337   Int_t isocount; //number of neighbours with 0 signal
338
339   Float_t d1[kDet][kMaxSMN][kMaxRow][kMaxCol];
340   
341   for(Int_t idet = 0; idet < kDet; idet++)
342     {
343       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
344         {
345           for(Int_t irow = 0; irow < kMaxRow; irow++)
346             {
347               for(Int_t icol = 0; icol < kMaxCol; icol++)
348                 {
349                   d1[idet][ismn][irow][icol] = 0.;
350                 }
351             }
352         }
353     }
354   
355   AliPMDRawStream rawStream(rawReader);
356
357   Int_t iddl = -1;
358   
359   Int_t numberofDDLs = 0;
360
361   while ((iddl = rawStream.DdlData(pmdddlcont)) >=0) {
362     numberofDDLs++;
363     
364     Int_t ientries = pmdddlcont->GetEntries();
365     
366     for (Int_t ient = 0; ient < ientries; ient++)
367       {
368         AliPMDddldata *pmdddl = (AliPMDddldata*)pmdddlcont->UncheckedAt(ient);
369         
370         Int_t idet = pmdddl->GetDetector();
371         Int_t ismn = pmdddl->GetSMN();
372         Int_t mcm = pmdddl->GetMCM();
373         //Int_t ichno = pmdddl->GetChannel();
374         Int_t irow = pmdddl->GetRow();
375         Int_t icol = pmdddl->GetColumn();
376         Int_t isig = pmdddl->GetSignal();
377         
378         // This is the protection not to crash the code 
379         
380         if(mcm == 0) continue;
381         if (irow < 0 || icol < 0 || irow > 47 || icol > 95) continue;
382         
383         // Pedestal subtraction
384         if(fHotFlag[idet][ismn][irow][icol] == 1.0) isig = 0;
385         
386         if (isig>0)
387           {
388             d1[idet][ismn][irow][icol] =
389               (Float_t) isig - fPedMeanRMS[idet][ismn][irow][icol];
390             
391             //printf("Signal_ped_subtracted=%f, pedestal=%f\n",d1[idet][ismn][irow][icol]),fPedMeanRMS[idet][ismn][irow][icol];
392             
393             fNhitCell[idet][ismn][irow][icol]++;     //cell hit frequency
394           }
395         
396
397       }
398     pmdddlcont->Delete();
399     }
400   
401   for(Int_t idet=0; idet < kDet; idet++)
402     {
403       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
404         {
405           for(Int_t irow = 0; irow < kMaxRow; irow++)
406             {
407               for(Int_t icol = 0; icol < kMaxCol; icol++)
408                 {
409                   if(d1[idet][ismn][irow][icol] > 0)
410                     {
411                       isocount = 0;
412                       for(Int_t ii = 0; ii < kCellNeighbour; ii++)
413                         {
414                           id1 = irow + neibx[ii];
415                           jd1 = icol + neiby[ii];
416                           if (id1 < 0) id1 = 0;
417                           if (id1 > kMaxRow-1) id1 = kMaxRow - 1;
418                           if (jd1 < 0) jd1 = 0;
419                           if (jd1 > kMaxCol-1) jd1 = kMaxCol - 1;
420                           if(d1[idet][ismn][id1][jd1] == 0)
421                             {
422                               isocount++;
423                               if(isocount == kCellNeighbour)
424                                 {
425                                   //fDetIso[idet] += d1[idet][ismn][irow][icol];
426                                   fSMIso[idet][ismn] += d1[idet][ismn][irow][icol];
427                                   fCellIso[idet][ismn][irow][icol] += d1[idet][ismn][irow][icol];
428                                   //fDetCount[idet]++;
429                                   fSMCount[idet][ismn]++;
430                                   fCellCount[idet][ismn][irow][icol]++;
431                                   
432                                 }
433                             }
434                         }  // neigh cell cond.
435                     }     // d>0 cond.
436                 }
437             }
438         }
439     }
440   
441   for(Int_t idet=0; idet < kDet; idet++)
442     {
443       for(Int_t ismn = 0; ismn < kMaxSMN; ismn++)
444         {
445           for(Int_t irow = 0; irow < kMaxRow; irow++)
446             {
447               for(Int_t icol = 0; icol < kMaxCol; icol++)
448                 {
449                   if(fNhitCell[idet][ismn][irow][icol]>0)
450                     {
451                       fCountSm[idet][ismn]  += 1;
452                       fTempnhit[idet][ismn] += fNhitCell[idet][ismn][irow][icol];
453                       fTempnhitSq[idet][ismn] += fNhitCell[idet][ismn][irow][icol]
454                         *fNhitCell[idet][ismn][irow][icol];
455                     }
456                 }
457             }
458         }
459     }
460   
461   if (numberofDDLs < kDDL) return kFALSE;
462   return kTRUE;
463   
464 }
465 // ------------------------------------------------------------------------ //
466 void AliPMDCalibGain::Analyse(TTree *gaintree, TTree *meantree)
467 {
468   // Calculates the mean
469   Int_t   det, sm, row, col;
470   Float_t gain;
471   Float_t cellmean = 0.;
472
473   Float_t modmean[2][24];
474
475   for (Int_t idet=0; idet < 2; idet++)
476     {
477       for (Int_t ism = 0; ism < 24; ism++)
478         {
479           modmean[idet][ism] = 0.;
480         }
481     }
482
483   gaintree->Branch("det",&det,"det/I");
484   gaintree->Branch("sm",&sm,"sm/I");
485   gaintree->Branch("row",&row,"row/I");
486   gaintree->Branch("col",&col,"col/I");
487   gaintree->Branch("gain",&gain,"gain/F");
488   
489   for(Int_t idet = 0; idet < kDet; idet++)
490     {
491       for(Int_t ism = 0; ism < kMaxSMN; ism++)
492         {
493           if (fSMCount[idet][ism] > 0)
494             modmean[idet][ism] = fSMIso[idet][ism]/fSMCount[idet][ism];
495           for(Int_t irow = 0; irow < kMaxRow; irow++)
496             {
497               for(Int_t icol = 0; icol < kMaxCol; icol++)
498                 {
499                   if (fCellCount[idet][ism][irow][icol] > 0.)
500                     {
501                       cellmean = fCellIso[idet][ism][irow][icol]/fCellCount[idet][ism][irow][icol];
502                     }
503                   det      = idet;
504                   sm       = ism;
505                   row      = irow;
506                   col      = icol;
507                   if (cellmean > 0.0 && fCellCount[idet][ism][irow][icol]>0.)
508                     {
509                       gain = cellmean/modmean[idet][ism];
510                     }
511                   else
512                     {
513                       gain = 0.;
514                     }
515                   //if(fCellCount[idet][ism][irow][icol]>0.) printf("CellCount =%f, gain= %f\n",fCellCount[idet][ism][irow][icol],gain);
516                   gaintree->Fill();
517                 }
518             }
519         }
520     }
521
522   Float_t smmean;
523
524   // Writing each module mean value
525   meantree->Branch("det",&det,"det/I");
526   meantree->Branch("sm",&sm,"sm/I");
527   meantree->Branch("smmean",&smmean,"row/F");
528   
529   for(Int_t idet = 0; idet < kDet; idet++)
530     {
531       for (Int_t ism = 0; ism < kMaxSMN; ism++)
532         {
533           det    = idet;
534           sm     = ism;
535           smmean = modmean[idet][ism];
536           meantree->Fill();
537         }
538     }
539
540 }
541 // ------------------------------------------------------------------------ //
542 void AliPMDCalibGain::FindHotCell(TTree *hottree, Float_t xvar)
543 {
544   // Calculates the mean
545   Int_t   det, sm, row, col;
546   Float_t flag;
547   Float_t meannhit;
548   Float_t meanSqnhit;
549   Float_t sigmanhit,nhitcut;
550
551   //Float_t xvar = 5.;
552
553   hottree->Branch("det",&det,"det/I");
554   hottree->Branch("sm",&sm,"sm/I");
555   hottree->Branch("row",&row,"row/I");
556   hottree->Branch("col",&col,"col/I");
557   hottree->Branch("flag",&flag,"flag/F");
558   
559   for(Int_t idet = 0; idet < kDet; idet++)
560     {
561       for(Int_t ism = 0; ism < kMaxSMN; ism++)
562         {
563           if (fCountSm[idet][ism]> 0)
564             {
565               meannhit   = fTempnhit[idet][ism]/fCountSm[idet][ism];
566               meanSqnhit = fTempnhitSq[idet][ism]/fCountSm[idet][ism];
567               sigmanhit  = sqrt(meanSqnhit-(meannhit*meannhit));
568               nhitcut    = meannhit + xvar*sigmanhit;
569
570               for(Int_t irow = 0; irow < kMaxRow; irow++)
571                 {
572                   for(Int_t icol = 0; icol < kMaxCol; icol++)
573                     {
574                       det      = idet;
575                       sm       = ism;
576                       row      = irow;
577                       col      = icol;
578                       
579                       if(fNhitCell[idet][ism][irow][icol] > nhitcut)
580                         {
581                           flag  = 1.0;
582                         }
583                       else
584                         {
585                           flag = 0.;
586                         }
587                       hottree->Fill();
588                     }
589                   
590                 }
591             }
592         }
593     }
594 }
595