Updated geometry including for the first time ACORDE
[u/mrichter/AliRoot.git] / EMCAL / AliCaloCalibPedestal.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15 /* $Id$ */
16
17 //________________________________________________________________________
18 //
19 // A help class for monitoring and calibration tools: MOOD, AMORE etc.,
20 // It can be created and used a la (ctor):
21 /*
22   //Create the object for making the histograms
23   fPedestals = new AliCaloCalibPedestal( fDetType );
24   // AliCaloCalibPedestal knows how many modules we have for PHOS or EMCAL
25   fNumModules = fPedestals->GetModules();
26 */
27 // fed an event:
28 //  fPedestals->ProcessEvent(fCaloRawStream);
29 // asked to draw histograms:
30 //  fPedestals->GetDeadMap(i)->Draw("col");
31 // or
32 //  fPedestals->GetPeakProfileHighGainRatio((i < fNumModules) ? i : fVisibleModule)->Draw("colz");
33 // etc.
34 // The pseudo-code examples above were from the first implementation in MOOD (summer 2007).
35 //________________________________________________________________________
36
37 #include "TH1.h"
38 #include "TFile.h"
39
40 #include <fstream>
41 #include <iostream>
42
43 #include "AliCaloRawStream.h"
44
45 //The include file
46 #include "AliCaloCalibPedestal.h"
47
48 ClassImp(AliCaloCalibPedestal)
49
50 using namespace std;
51
52 // ctor; initialize everything in order to avoid compiler warnings
53 AliCaloCalibPedestal::AliCaloCalibPedestal(kDetType detectorType) :  
54   TObject(),
55   fPedestalLowGain(),
56   fPedestalHighGain(),
57   fPeakMinusPedLowGain(),
58   fPeakMinusPedHighGain(),
59   fPedestalLowGainDiff(),
60   fPedestalHighGainDiff(),
61   fPeakMinusPedLowGainDiff(),
62   fPeakMinusPedHighGainDiff(),
63   fPedestalLowGainRatio(),
64   fPedestalHighGainRatio(),
65   fPeakMinusPedLowGainRatio(),
66   fPeakMinusPedHighGainRatio(),
67   fDeadMap(),
68   fDeadTowers(0),
69   fNewDeadTowers(0),
70   fResurrectedTowers(0),
71   fReference(0),
72   fDetType(kNone),
73   fColumns(0),
74   fRows(0),
75   fModules(0),
76   fRunNumber(-1)
77 {
78   //Default constructor. First we set the detector-type related constants.
79   if (detectorType == kPhos) {
80     fColumns = fgkPhosCols;
81     fRows = fgkPhosRows;
82     fModules = fgkPhosModules;
83   } 
84   else {
85     //We'll just trust the enum to keep everything in line, so that if detectorType
86     //isn't kPhos then it is kEmCal. Note, however, that this is not necessarily the
87     //case, if someone intentionally gives another number
88     fColumns = fgkEmCalCols;
89     fRows = fgkEmCalRows;
90     fModules = fgkEmCalModules;
91   } 
92   fDetType = detectorType;
93  
94   //Then, loop for the requested number of modules
95   TString title, name;
96   for (int i = 0; i < fModules; i++) {
97     //Pedestals, low gain
98     name = "hPedlowgain";
99     name += i;
100     title = "Pedestals, low gain, module ";
101     title += i; 
102     fPedestalLowGain.Add(new TProfile2D(name, title,
103                                         fColumns, 0.0, fColumns, 
104                                         fRows, -fRows, 0.0));
105   
106     //Pedestals, high gain
107     name = "hPedhighgain";
108     name += i;
109     title = "Pedestals, high gain, module ";
110     title += i; 
111     fPedestalHighGain.Add(new TProfile2D(name, title,
112                                          fColumns, 0.0, fColumns, 
113                                          fRows, -fRows, 0.0));
114   
115     //Peak-Pedestals, low gain
116     name = "hPeakMinusPedlowgain";
117     name += i;
118     title = "Peak-Pedestal, low gain, module ";
119     title += i; 
120     fPeakMinusPedLowGain.Add(new TProfile2D(name, title,
121                                             fColumns, 0.0, fColumns, 
122                                             fRows, -fRows, 0.0));
123   
124     //Peak-Pedestals, high gain
125     name = "hPeakMinusPedhighgain";
126     name += i;
127     title = "Peak-Pedestal, high gain, module ";
128     title += i; 
129     fPeakMinusPedHighGain.Add(new TProfile2D(name, title,
130                                              fColumns, 0.0, fColumns, 
131                                              fRows, -fRows, 0.0));
132   
133     name = "hDeadMap";
134     name += i;
135     title = "Dead map, module ";
136     title += i;
137     fDeadMap.Add(new TH2D(name, title, fColumns, 0.0, fColumns, 
138                           fRows, -fRows, 0.0));
139   
140   }//end for nModules create the histograms
141  
142   //Compress the arrays, in order to remove the empty objects (a 16 slot array is created by default)
143   fPedestalLowGain.Compress();
144   fPedestalHighGain.Compress();
145   fPeakMinusPedLowGain.Compress();
146   fPeakMinusPedHighGain.Compress();
147   fDeadMap.Compress();
148   //Make them the owners of the profiles, so we don't need to care about deleting them
149   //fPedestalLowGain.SetOwner();
150   //fPedestalHighGain.SetOwner();
151   //fPeakMinusPedLowGain.SetOwner();
152   //fPeakMinusPedHighGain.SetOwner();
153   
154 }
155
156 // dtor
157 //_____________________________________________________________________
158 AliCaloCalibPedestal::~AliCaloCalibPedestal()
159 {
160   if (fReference) delete fReference;//Delete the reference object, if it has been loaded
161   //TObjArray will delete the histos/profiles when it is deleted.
162 }
163
164 // copy ctor
165 //_____________________________________________________________________
166 AliCaloCalibPedestal::AliCaloCalibPedestal(const AliCaloCalibPedestal &ped) :
167   TObject(ped),
168   fPedestalLowGain(),
169   fPedestalHighGain(),
170   fPeakMinusPedLowGain(),
171   fPeakMinusPedHighGain(),
172   fPedestalLowGainDiff(),
173   fPedestalHighGainDiff(),
174   fPeakMinusPedLowGainDiff(),
175   fPeakMinusPedHighGainDiff(),
176   fPedestalLowGainRatio(),
177   fPedestalHighGainRatio(),
178   fPeakMinusPedLowGainRatio(),
179   fPeakMinusPedHighGainRatio(),
180   fDeadMap(),
181   fDeadTowers(ped.GetDeadTowerCount()),
182   fNewDeadTowers(ped.GetDeadTowerNew()),
183   fResurrectedTowers(ped.GetDeadTowerResurrected()),
184   fReference( 0 ), //! note that we do not try to copy the reference info here
185   fDetType(ped.GetDetectorType()),
186   fColumns(ped.GetColumns()),
187   fRows(ped.GetRows()),
188   fModules(ped.GetModules()),
189   fRunNumber(ped.GetRunNumber())
190 {
191   // Then the ObjArray ones; we add the histograms rather than trying TObjArray = assignment
192   //DS: this has not really been tested yet..
193   for (int i = 0; i < fModules; i++) {
194     fPedestalLowGain.Add( ped.GetPedProfileLowGain(i) );
195     fPedestalHighGain.Add( ped.GetPedProfileHighGain(i) );
196     fPeakMinusPedLowGain.Add( ped.GetPeakProfileLowGain(i) );
197     fPeakMinusPedHighGain.Add( ped.GetPeakProfileHighGain(i) );
198
199     fDeadMap.Add( ped.GetDeadMap(i) );  
200   }//end for nModules 
201  
202   //Compress the arrays, in order to remove the empty objects (a 16 slot array is created by default)
203   fPedestalLowGain.Compress();
204   fPedestalHighGain.Compress();
205   fPeakMinusPedLowGain.Compress();
206   fPeakMinusPedHighGain.Compress();
207   fDeadMap.Compress();
208 }
209
210 // assignment operator; use copy ctor to make life easy..
211 //_____________________________________________________________________
212 AliCaloCalibPedestal& AliCaloCalibPedestal::operator = (const AliCaloCalibPedestal &source)
213 {
214   // assignment operator; use copy ctor
215   if (&source == this) return *this;
216
217   new (this) AliCaloCalibPedestal(source);
218   return *this;
219 }
220
221 //_____________________________________________________________________
222 void AliCaloCalibPedestal::Reset()
223 {
224   // Reset all arrays/histograms
225   for (int i = 0; i < fModules; i++) {
226     GetPedProfileLowGain(i)->Reset();
227     GetPedProfileHighGain(i)->Reset();
228     GetPeakProfileLowGain(i)->Reset();
229     GetPeakProfileHighGain(i)->Reset();
230     GetDeadMap(i)->Reset();
231     
232     if (!fPedestalLowGainDiff.IsEmpty()) {
233       //This means that the comparison profiles have been created.
234   
235       GetPedProfileLowGainDiff(i)->Reset();
236       GetPedProfileHighGainDiff(i)->Reset();
237       GetPeakProfileLowGainDiff(i)->Reset();
238       GetPeakProfileHighGainDiff(i)->Reset();
239       
240       GetPedProfileLowGainRatio(i)->Reset();
241       GetPedProfileHighGainRatio(i)->Reset();
242       GetPeakProfileLowGainRatio(i)->Reset();
243       GetPeakProfileHighGainRatio(i)->Reset();
244     }
245   }
246   fDeadTowers = 0;
247   fNewDeadTowers = 0;
248   fResurrectedTowers = 0;
249  
250   //To think about: should fReference be deleted too?... let's not do it this time, at least...
251 }
252
253 //_____________________________________________________________________
254 Bool_t AliCaloCalibPedestal::ProcessEvent(AliCaloRawStream *in)
255
256   // Method to process=analyze one event in the data stream
257   if (!in) return kFALSE; //Return right away if there's a null pointer
258   
259   int sample, i = 0; //The sample temp, and the sample number in current event.
260   int max = fgkSampleMin, min = fgkSampleMax;//Use these for picking the pedestal
261   int gain = 0;
262   
263   while (in->Next()) {
264     sample = in->GetSignal(); //Get the adc signal
265     if (sample < min) min = sample;
266     if (sample > max) max = sample;
267     i++;
268     if ( i >= in->GetTimeLength()) {
269       //If we're here then we're done with this tower
270       gain = 1 - in->IsLowGain();
271       
272       int arrayPos = in->GetModule(); //The modules are numbered starting from 0
273       if (arrayPos >= fModules) {
274         //TODO: return an error message, if appopriate (perhaps if debug>0?)
275         return kFALSE;
276       } 
277     
278       //Debug
279       if (arrayPos < 0 || arrayPos >= fModules) {
280         printf("Oh no: arrayPos = %i.\n", arrayPos); 
281       }
282
283       //NOTE: coordinates are (column, row) for the profiles
284       if (gain == 0) {
285         //fill the low gain histograms
286         ((TProfile2D*)fPedestalLowGain[arrayPos])->Fill(in->GetColumn(), -in->GetRow() - 1, min);
287         ((TProfile2D*)fPeakMinusPedLowGain[arrayPos])->Fill(in->GetColumn(), -in->GetRow() - 1, max - min);
288       } 
289       else {//fill the high gain ones
290         ((TProfile2D*)fPedestalHighGain[arrayPos])->Fill(in->GetColumn(), -in->GetRow() - 1, min);
291         ((TProfile2D*)fPeakMinusPedHighGain[arrayPos])->Fill(in->GetColumn(), -in->GetRow() - 1, max - min);
292       }//end if gain
293       
294       max = fgkSampleMin; min = fgkSampleMax;
295       i = 0;
296     
297     }//End if end of tower
298    
299   }//end while, of stream
300   
301   return kTRUE;
302 }
303
304 //_____________________________________________________________________
305 Bool_t AliCaloCalibPedestal::SaveHistograms(TString fileName, Bool_t saveEmptyHistos)
306 {
307   //Saves all the histograms (or profiles, to be accurate) to the designated file
308   
309   TFile destFile(fileName, "recreate");
310   
311   if (destFile.IsZombie()) {
312     return kFALSE;
313   }
314   
315   destFile.cd();
316   
317   for (int i = 0; i < fModules; i++) {
318     if( ((TProfile2D *)fPeakMinusPedLowGain[i])->GetEntries() || saveEmptyHistos) {
319       fPeakMinusPedLowGain[i]->Write();
320     }
321     if( ((TProfile2D *)fPeakMinusPedHighGain[i])->GetEntries() || saveEmptyHistos) { 
322       fPeakMinusPedHighGain[i]->Write();
323     }
324     if( ((TProfile2D *)fPedestalLowGain[i])->GetEntries() || saveEmptyHistos) {
325       fPedestalLowGain[i]->Write();
326     }
327     if( ((TProfile2D *)fPedestalHighGain[i])->GetEntries() || saveEmptyHistos) {
328       fPedestalHighGain[i]->Write();
329     }
330   } 
331   
332   destFile.Close();
333   
334   return kTRUE;
335 }
336
337 //_____________________________________________________________________
338 Bool_t AliCaloCalibPedestal::LoadReferenceCalib(TString fileName, TString objectName)
339 {
340   
341   //Make sure that the histograms created when loading the object are not destroyed as the file object is destroyed
342   TH1::AddDirectory(kFALSE);
343   
344   TFile *sourceFile = new TFile(fileName);
345   if (sourceFile->IsZombie()) {
346     return kFALSE;//We couldn't load the reference
347   }
348
349   if (fReference) delete fReference;//Delete the reference object, if it already exists
350   fReference = 0;
351   
352   fReference = (AliCaloCalibPedestal*)sourceFile->Get(objectName);
353  
354   if (!fReference || !(fReference->InheritsFrom(AliCaloCalibPedestal::Class())) || (fReference->GetDetectorType() != fDetType)) {
355     if (fReference) delete fReference;//Delete the object, in case we had an object of the wrong type
356     fReference = 0;
357     return kFALSE;
358   }
359         
360   delete sourceFile;
361  
362   //Reset the histogram ownership behaviour. NOTE: a better workaround would be good, since this may accidentally set AddDirectory to true, even
363   //if we are called by someone who has set it to false...
364   TH1::AddDirectory(kTRUE);
365  
366   return kTRUE;//We succesfully loaded the object
367 }
368
369 //_____________________________________________________________________
370 void AliCaloCalibPedestal::ValidateComparisonProfiles()
371 {
372   //Make sure the comparison histos exist
373   if (!fPedestalLowGainDiff.IsEmpty()) return; //The profiles already exist. We just check one, because they're all created at
374   //the same time
375                                                 
376                                                 
377   //Then, loop for the requested number of modules
378   TString title, name;
379   for (int i = 0; i < fModules; i++) {
380     //Pedestals, low gain
381     name = "hPedlowgainDiff";
382     name += i;
383     title = "Pedestals difference, low gain, module ";
384     title += i; 
385     fPedestalLowGainDiff.Add(new TProfile2D(name, title,
386                                             fColumns, 0.0, fColumns, 
387                                             fRows, -fRows, 0.0));
388   
389     //Pedestals, high gain
390     name = "hPedhighgainDiff";
391     name += i;
392     title = "Pedestals difference, high gain, module ";
393     title += i; 
394     fPedestalHighGainDiff.Add(new TProfile2D(name, title,
395                                              fColumns, 0.0, fColumns, 
396                                              fRows, -fRows, 0.0));
397   
398     //Peak-Pedestals, low gain
399     name = "hPeakMinusPedlowgainDiff";
400     name += i;
401     title = "Peak-Pedestal difference, low gain, module ";
402     title += i; 
403     fPeakMinusPedLowGainDiff.Add(new TProfile2D(name, title,
404                                                 fColumns, 0.0, fColumns, 
405                                                 fRows, -fRows, 0.0));
406   
407     //Peak-Pedestals, high gain
408     name = "hPeakMinusPedhighgainDiff";
409     name += i;
410     title = "Peak-Pedestal difference, high gain, module ";
411     title += i; 
412     fPeakMinusPedHighGainDiff.Add(new TProfile2D(name, title,
413                                                  fColumns, 0.0, fColumns, 
414                                                  fRows, -fRows, 0.0));
415   
416     //Pedestals, low gain
417     name = "hPedlowgainRatio";
418     name += i;
419     title = "Pedestals ratio, low gain, module ";
420     title += i; 
421     fPedestalLowGainRatio.Add(new TProfile2D(name, title,
422                                              fColumns, 0.0, fColumns, 
423                                              fRows, -fRows, 0.0));
424   
425     //Pedestals, high gain
426     name = "hPedhighgainRatio";
427     name += i;
428     title = "Pedestals ratio, high gain, module ";
429     title += i; 
430     fPedestalHighGainRatio.Add(new TProfile2D(name, title,
431                                               fColumns, 0.0, fColumns, 
432                                               fRows, -fRows, 0.0));
433   
434     //Peak-Pedestals, low gain
435     name = "hPeakMinusPedlowgainRatio";
436     name += i;
437     title = "Peak-Pedestal ratio, low gain, module ";
438     title += i; 
439     fPeakMinusPedLowGainRatio.Add(new TProfile2D(name, title,
440                                                  fColumns, 0.0, fColumns, 
441                                                  fRows, -fRows, 0.0));
442   
443     //Peak-Pedestals, high gain
444     name = "hPeakMinusPedhighgainRatio";
445     name += i;
446     title = "Peak-Pedestal ratio, high gain, module ";
447     title += i; 
448     fPeakMinusPedHighGainRatio.Add(new TProfile2D(name, title,
449                                                   fColumns, 0.0, fColumns, 
450                                                   fRows, -fRows, 0.0));
451     
452   }//end for nModules create the histograms
453 }
454
455 //_____________________________________________________________________
456 void AliCaloCalibPedestal::ComputeDiffAndRatio()
457 {
458   // calculate differences and ratios relative to a reference
459   ValidateComparisonProfiles();//Make sure the comparison histos exist
460  
461   if (!fReference) {
462     return;//Return if the reference object isn't loaded
463   }
464
465   for (int i = 0; i < fModules; i++) {
466     //Compute the ratio of the histograms
467     
468     ((TProfile2D*)fPedestalLowGainRatio[i])->Divide(GetPedProfileLowGain(i), fReference->GetPedProfileLowGain(i), 1.0, 1.0);
469     ((TProfile2D*)fPedestalHighGainRatio[i])->Divide(GetPedProfileHighGain(i), fReference->GetPedProfileHighGain(i), 1.0, 1.0);
470     ((TProfile2D*)fPeakMinusPedLowGainRatio[i])->Divide(GetPeakProfileLowGain(i), fReference->GetPeakProfileLowGain(i), 1.0, 1.0);
471     ((TProfile2D*)fPeakMinusPedHighGainRatio[i])->Divide(GetPeakProfileHighGain(i), fReference->GetPeakProfileHighGain(i), 1.0, 1.0);
472   
473     //For computing the difference, we cannot simply do TProfile2D->Add(), because that subtracts the sum of all entries,
474     //which means that the mean of the new profile will not be the difference of the means. So do it by hand:
475     for (int j = 0; j <= fColumns; j++) {
476       for (int k = 0; k <= fRows; k++) {
477         int bin = ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->GetBin(j+1, k+1);//Note that we assume here that all histos have the same structure...
478         double diff = fReference->GetPeakProfileHighGain(i)->GetBinContent(bin) - GetPeakProfileHighGain(i)->GetBinContent(bin);
479         ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinContent(j+1, k+1, diff);
480         ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinEntries(bin, 1);
481
482         diff = fReference->GetPeakProfileLowGain(i)->GetBinContent(bin) - GetPeakProfileLowGain(i)->GetBinContent(bin);
483         ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinContent(j+1, k+1, diff);
484         ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinEntries(bin, 1);
485     
486         diff = fReference->GetPedProfileHighGain(i)->GetBinContent(bin) - GetPedProfileHighGain(i)->GetBinContent(bin);
487         ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinContent(j+1, k+1, diff);
488         ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinEntries(bin, 1);
489
490         diff = fReference->GetPedProfileLowGain(i)->GetBinContent(bin) - GetPedProfileLowGain(i)->GetBinContent(bin);
491         ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinContent(j+1, k+1, diff);
492         ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinEntries(bin, 1);
493        
494       } // rows
495     } // columns
496     
497   } // modules
498  
499 }
500
501 //_____________________________________________________________________
502 void AliCaloCalibPedestal::ComputeDeadTowers(int threshold, const char * deadMapFile)
503 {
504   //Computes the number of dead towers etc etc into memory, after this you can call the GetDead... -functions
505   int countTot = 0;
506   int countNew = 0;
507   int countRes = 0;
508   ofstream * fout = 0;
509   ofstream * diff = 0;
510   char name[512];//Quite a long temp buffer, just in case the filename includes a path
511   
512   if (deadMapFile) {
513     snprintf(name, 512, "%s.txt", deadMapFile);
514     fout = new ofstream(name);
515     snprintf(name, 512, "%sdiff.txt", deadMapFile);
516     diff = new ofstream(name);
517     if (!fout->is_open()) {
518       delete fout;
519       fout = 0;//Set the pointer to empty if the file was not opened
520     }
521     if (!diff->is_open()) {
522       delete diff;
523       fout = 0;//Set the pointer to empty if the file was not opened
524     }
525   }
526  
527   for (int i = 0; i < fModules; i++) {
528     if (GetPeakProfileHighGain(i)->GetEntries() > 0) { //don't care about empty histos
529       for (int j = 1; j <= fColumns; j++) {
530         for (int k = 1; k <= fRows; k++) {
531
532           if (GetPeakProfileHighGain(i)->GetBinContent(j, k) < threshold) {//It's dead
533             countTot++;//One more dead total
534             if (fout) {
535               (*fout) << i << " " 
536                       << (fRows - k) << " " 
537                       << (j-1) << " " 
538                       << "1" << " " 
539                       << "0" << endl;//Write the status to the deadmap file, if the file is open.
540             }
541             
542             if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) >= threshold) {
543               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kRecentlyDeceased); 
544               countNew++;//This tower wasn't dead before!
545               if (diff) {
546                 ( *diff) << i << " " 
547                          << (fRows - k) << " " 
548                          << (j - 1) << " " 
549                          << "1" << " " 
550                          << "0" << endl;//Write the status to the deadmap difference file, if the file is open.
551               }
552             } 
553             else {
554               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kDead);//This has been dead before. Nothing new         
555             }
556           } 
557           else { //It's ALIVE!!
558             //Don't bother with writing the live ones.
559             //if (fout)
560             //  (*fout) << i << " " 
561             //     << (fRows - k) << " " 
562             //     << (j - 1) << " " 
563             //     << "1" << " " 
564             //     << "1" << endl;//Write the status to the deadmap file, if the file is open.
565             if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) < threshold) {
566               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kResurrected);
567               countRes++; //This tower was dead before => it's a miracle! :P
568               if (diff) {
569                 (*diff) << i << " " 
570                         << (fRows - k) << " " 
571                         << (j - 1) << " " 
572                         << "1" << " " 
573                         << "1" << endl;//Write the status to the deadmap difference file, if the file is open.
574               }
575             } 
576             else {
577               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kAlive);
578             }
579           }
580             
581         }//end for k/rows
582       }//end for j/columns
583     }//end if GetEntries >= 0
584   
585   }//end for modules
586  
587  if (fout) {
588    fout->close();
589    delete fout;
590  }
591  
592  fDeadTowers = countTot;
593  fNewDeadTowers = countNew;
594  fResurrectedTowers = countRes;
595 }
596