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