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