Preprocessor classes for SPD (Paul Nilsson)
[u/mrichter/AliRoot.git] / ITS / AliITSPreprocessorSPD.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
16 /*
17 $Log$
18 */
19
20 ///////////////////////////////////////////////////////////////////////////
21 // AliITSPreprocessorSPD implementation by P. Nilsson 2005
22 // MAIN AUTHOR/CONTACT: Paul.Nilsson@cern.ch
23 //
24 // The purpose of this class is to provide algorithms for identification
25 // of "bad channels" such as dead channels and noisy channels in the SPD
26 //
27 // Examples on how to use this class can be found in the macros:
28 //
29 // .findNoisyChannels.C (Locate and store noisy channels in the CDB)
30 // .readNoisyChannels.C (Read noisy channels from the CDB)
31 // .findDeadChannels.C  (Locate and store dead channels in the CDB)
32 // .readDeadChannels.C  (Read dead channels from the CDB)
33 ///////////////////////////////////////////////////////////////////////////
34
35 #include "AliITSPreprocessorSPD.h"
36
37 ClassImp(AliITSPreprocessorSPD)
38
39
40 //__________________________________________________________________________
41 AliITSPreprocessorSPD::AliITSPreprocessorSPD(void):
42   fITSLoader(0x0),
43   fRunLoader(0x0),
44   fThresholdRatio(5.),
45   fThreshold(5),
46   fMaximumNumberOfEvents(1000000),
47   fHighestModuleNumber(0),
48   fSelectedAlgorithm(kOptimizedForRealData),
49   fGeometryMode(kALICEGeometry),
50   fNumberOfBadChannels(0),
51   fInit(kFALSE),
52   fVMEMode(kFALSE),
53   fDigitsHistogram(0),
54   fBadChannelsObjArray(0),
55   fBadChannelsIntArray(0),
56   fBadChannelsIndexArray(0),
57   fBadChannelsContainer(0)
58 {
59   // Default constructor for the SPD preprocessor
60   //
61   // Initialization has to be done by hand using Set* methods and the Open method
62   //
63   // Input :
64   // Output:
65   // Return: <empty/uninitialized AliITSPreprocessorSPD object>
66 }
67
68
69 //__________________________________________________________________________
70 AliITSPreprocessorSPD::AliITSPreprocessorSPD(const char *fileName, const char *mode,
71                                              const char *fileNameg, const Int_t maxNumberOfEvents):
72   fITSLoader(0x0),
73   fRunLoader(0x0),
74   fInit(kFALSE)
75 {
76   // Standard constructor for the SPD preprocessor
77   //
78   // Input : Name of digit file
79   // Output: fInit
80   // Return: Initialized SPD preprocessor object
81
82   // Initialize
83   AliITSPreprocessorSPD::Init();
84   AliITSPreprocessorSPD::SetMaximumNumberOfEvents(maxNumberOfEvents);
85
86   // Open and read the galice and digit files
87   if (!AliITSPreprocessorSPD::Open(fileName, mode, fileNameg))
88     {
89       AliError("Failed to open file");
90       fInit = kFALSE;
91     }
92 }
93
94
95 //__________________________________________________________________________
96 AliITSPreprocessorSPD::AliITSPreprocessorSPD(const AliITSPreprocessorSPD &prep) :
97   TTask(prep)
98 {
99   // Default copy constructor
100   // Notice that only pointer addresses are copied!
101   // Memory allocations of new objects are not done.
102
103   fITSLoader = prep.fITSLoader;
104   fRunLoader = prep.fRunLoader;
105   fThresholdRatio = prep.fThresholdRatio;
106   fThreshold = prep.fThreshold;
107   fMaximumNumberOfEvents = prep.fMaximumNumberOfEvents;
108   fHighestModuleNumber = prep.fHighestModuleNumber;
109   fSelectedAlgorithm = prep.fSelectedAlgorithm;
110   fGeometryMode = prep.fGeometryMode;
111   fNumberOfBadChannels = prep.fNumberOfBadChannels;
112   fInit = prep.fInit;
113   fVMEMode = prep.fVMEMode;
114   fDigitsHistogram = prep.fDigitsHistogram;
115   fBadChannelsObjArray = prep.fBadChannelsObjArray;
116   fBadChannelsIntArray = prep.fBadChannelsIntArray;
117   fBadChannelsIndexArray = prep.fBadChannelsIndexArray;
118   fBadChannelsContainer = prep.fBadChannelsContainer;
119 }
120
121
122 //__________________________________________________________________________
123 AliITSPreprocessorSPD& AliITSPreprocessorSPD::operator=(const AliITSPreprocessorSPD &prep)
124 {
125   // Assignment operator
126
127   AliError("Not implemented");
128
129   if (this != &prep) { } // Do not delete this line
130
131   return *this;
132 }
133
134
135 //__________________________________________________________________________
136 AliITSPreprocessorSPD::~AliITSPreprocessorSPD(void)
137 {
138   // Destructor for the SPD preprocessor
139
140   if (fDigitsHistogram)
141     {
142       // try fDigitsHistogram->Delete(); if the following crashes
143       for (UInt_t module = 0; module < fNumberOfModules; module++)
144         {
145           (*fDigitsHistogram)[module]->Delete();
146         }
147       delete fDigitsHistogram;
148       fDigitsHistogram = 0;
149     }
150
151   if (fNumberOfBadChannels)
152     {
153       delete [] fNumberOfBadChannels;
154       fNumberOfBadChannels = 0;
155     }
156
157   if (fBadChannelsIntArray)
158     {
159       delete [] fBadChannelsIntArray;
160       fBadChannelsIntArray = 0;
161     }
162
163   if (fBadChannelsIndexArray)
164     {
165       delete [] fBadChannelsIndexArray;
166       fBadChannelsIndexArray = 0;
167     }
168
169   if (fBadChannelsObjArray)
170     {
171       fBadChannelsObjArray->Delete();
172       delete fBadChannelsObjArray;
173       fBadChannelsObjArray = 0;
174     }
175
176   delete fRunLoader;
177   fRunLoader = 0;
178
179   delete fBadChannelsContainer;
180   fBadChannelsContainer = 0;
181 }
182
183
184 //__________________________________________________________________________
185 void AliITSPreprocessorSPD::Init(void)
186 {
187   // Initialization of the SPD preprocessor
188   //
189   // Input : (void)
190   // Output: Several logistical variables
191   // Return: (void)
192
193   // Default maximum number of events per histogram
194   fMaximumNumberOfEvents = 1000000;
195
196   // Default noisy channel removal algorithm
197   fSelectedAlgorithm = kOptimizedForRealData;
198
199   // Default noisy channel threshold and threshold ratio
200   // (threshold for current bin content divided by the average neighboring bin contents)
201   fThreshold = 10;
202   fThresholdRatio = 5.;
203
204   // Set default geometry mode (ALICE geometry). This also sets fNumberOfModules
205   AliITSPreprocessorSPD::SetGeometryMode(kALICEGeometry);
206
207   // The highest module number with found bad channels
208   fHighestModuleNumber = 0;
209
210   // Assume input is not from a VME file
211   fVMEMode = kFALSE;
212
213   // Initialization is complete
214   fInit = kTRUE;
215 }
216
217
218 //__________________________________________________________________________
219 void AliITSPreprocessorSPD::SetGeometryMode(UInt_t mode)
220 {
221   // Set the geometry mode
222   //
223   // Input : Geometry mode (either kTestBeamGeometry or kALICEGeometry)
224   // Output: 
225   // Return: (void)
226
227   fGeometryMode = mode;
228
229   // In case of an input VME file, the number of modules has already been fixed.
230   // Do not try to change it
231   if (!fVMEMode)
232     {
233       if (mode == kALICEGeometry)
234         {
235           fNumberOfModules = kNumberOfSPDModules;
236         }
237       else if (mode == kTestBeamGeometry)
238         {
239           fNumberOfModules = kNumberOfTestBeamSPDModules;
240         }
241       else
242         {
243           AliError("Unknown geometry mode, defaults to ALICE geometry");
244           fNumberOfModules = kNumberOfSPDModules;
245         }
246     }
247 }
248
249
250 //__________________________________________________________________________
251 void AliITSPreprocessorSPD::SetFakeNoisyChannel(Int_t module, Int_t column, Int_t row)
252 {
253   // Introduce a fake noisy channel in the hit histograms
254   //
255   // Input : Module, column and row numbers
256   // Output: Updated hit histograms
257   // Return: (void)
258
259   if ((UInt_t)module < fNumberOfModules)
260     {
261       ((TH2F*)(*fDigitsHistogram)[module])->Fill(column, row, 1000000);
262     }
263   else
264     {
265       AliError("No such module number");
266     }
267 }
268
269
270 //__________________________________________________________________________
271 Bool_t AliITSPreprocessorSPD::Open(const char *fileName, const char *mode, const char *fileNameg)
272 {
273   // Open digit file
274   //
275   // Input : Name and mode of ITS digit file, name of galice file
276   // Output: Digits
277   // Return: kTRUE if loaders succeed
278
279   Bool_t status = kFALSE;
280   Bool_t status0 = kFALSE;
281   Bool_t status1 = kFALSE;
282   Bool_t status2 = kFALSE;
283   Bool_t status3 = kFALSE;
284
285   // Only proceed if initialization has been done
286   if (fInit)
287     {
288       TString m = (TString)mode;
289       if (m == "daq" || m == "DAQ")
290         {
291           // Open the data file and get the run loader
292           fRunLoader = AliRunLoader::Open(fileNameg);
293           if (fRunLoader)
294             {
295               // Get the gAlice object
296               status0 = AliITSPreprocessorSPD::GetgAlice();
297
298               // Get the ITS digits
299               if (status0) status1 = AliITSPreprocessorSPD::GetITSDigits(fileName);
300
301               // Create the test beam object
302               if (status1) status2 = AliITSPreprocessorSPD::CreateGeometryObj();
303
304               // Fill histograms with DAQ digits
305               if (status2) status3 = AliITSPreprocessorSPD::FillHistograms();
306
307               status = status0 & status1 & status2 & status3;
308             }
309           else
310             {
311               AliError("Failed to get the run loader");
312             }
313         }
314       else if (m == "vme" || m == "VME")
315         {
316           // Open the VME file that contains histograms with fired channels as read by the stand-alone VME system
317           TFile *vmeFile = TFile::Open(fileName);
318
319           if (!vmeFile->IsOpen())
320             {
321               AliError("Could not open VME input file");
322             }
323           else
324             {
325               // Get the histograms from the VME file that contains all fired channels
326               status0 = AliITSPreprocessorSPD::GetVMEHistograms(vmeFile);
327
328               // Create the test beam object
329               if (status0) status1 = AliITSPreprocessorSPD::CreateGeometryObj();
330             }
331
332           // This boolean will be used to override any attempts of changing the number of modules by the user
333           // with the SetGeometryMode method. For VME files the number of modules will entirely be determined by
334           // the input VME root file, i.e. by the number of histograms in this file
335           fVMEMode = kTRUE;
336
337           status = status0 & status1;
338         }
339       else
340         {
341           AliError("Unknown filetype - assuming DAQ file");
342         }
343
344       // At this stage, the final number of modules will be known (safe to define arrays)
345       // In case data is read from a VME root file, it will not be known before
346       if (status)
347         {
348           // Create the arrays for bookkeeping and storing the noisy channels
349           if (!fBadChannelsObjArray)
350             {
351               fBadChannelsObjArray = new TObjArray();
352             }
353           if (!fBadChannelsIndexArray)
354             {
355               // This array will contain the begin and end indices for each module, i.e. where to begin
356               // and stop reading the fBadChannelsObjArray for a certain module.
357               // The last position of the array will contain the end index of the last module
358               fBadChannelsIndexArray = new Int_t[fNumberOfModules + 1];
359               for (UInt_t module = 0; module < fNumberOfModules + 1; module++)
360                 {
361                   fBadChannelsIndexArray[module] = 0;
362                 }
363             }
364         }
365     }
366   else
367     {
368       AliError("SPD preprocessor not initialized. Can't load digits");
369     }
370
371   return status;
372 }
373
374
375 //__________________________________________________________________________
376 Bool_t AliITSPreprocessorSPD::GetgAlice(void)
377 {
378   // Get the gAlice object
379   //
380   // Input : (void)
381   // Output: gAlice object
382   // Return: kTRUE if the gAlice object was found
383
384   Bool_t status = kTRUE;
385
386   // Load gAlice
387   fRunLoader->LoadgAlice();
388   if (!fRunLoader->GetAliRun())
389     {
390       AliError("gAlice not found on file. Aborting.");
391       status = kFALSE;
392     }
393
394   return status;
395 }
396
397
398 //__________________________________________________________________________
399 Bool_t AliITSPreprocessorSPD::GetVMEHistograms(TFile *vmeFile)
400 {
401   // Get pre-filled digit histograms from input VME file
402   //
403   // Input : pointer to VME file
404   // Output: 
405   // Return: kTRUE if histograms are found on file
406
407   Bool_t status = kFALSE;
408
409   // Get the file directory
410   TDirectory *dir = (TDirectory *)vmeFile;
411
412   // Get the number of keys (histograms in this case corresponding to modules/ladders)
413   fNumberOfModules = dir->GetNkeys();
414   if ((fNumberOfModules > 0) && (fNumberOfModules <= kNumberOfSPDModules))
415     {
416       status = kTRUE;
417
418       // Create bad channel histograms
419       fDigitsHistogram = new TObjArray(fNumberOfModules);
420
421       // Create a key iterator
422       TIter nextkey(dir->GetListOfKeys());
423       TKey *key = 0;
424
425       // Loop over all objects, read them in to memory one by one
426       UInt_t module = 0;
427       while ( (key = (TKey *)nextkey()) )
428         {
429           (*fDigitsHistogram)[module++] = (TH2F *)key->ReadObj();
430
431           // For safety
432           if (module > fNumberOfModules)
433             {
434               AliError("Wrong number of keys in VME file");
435               status = kFALSE;
436               break;
437             }
438         }
439     }
440   else
441     {
442       AliError("Wrong number of histograms in VME file");
443     }
444   return status;
445 }
446
447 //__________________________________________________________________________
448 Bool_t AliITSPreprocessorSPD::GetITSDigits(const char *fileName)
449 {
450   // Get the ITS digits
451   //
452   // Input : name of digits file
453   // Output: fITSLoader object, ITS digits
454   // Return: kTRUE if loader succeed
455
456   Bool_t status = kTRUE;
457
458   // Load the gAlice and the header
459   fRunLoader->LoadgAlice();
460   fRunLoader->LoadHeader();
461
462   // Get the ITS loader
463   fITSLoader = (AliITSLoader*) fRunLoader->GetLoader("ITSLoader");
464   if (!fITSLoader)
465     {
466       AliError("ITS loader not found");
467       status = kFALSE;
468     }
469   else
470     {
471       // Open the digits file
472       fITSLoader->SetDigitsFileName(fileName);
473     }
474
475   return status;
476 }
477
478
479 //__________________________________________________________________________
480 TClonesArray* AliITSPreprocessorSPD::CreateDigitsArray(void) const
481 {
482   // Creation of the digits array
483   //
484   // Input : (void)
485   // Output:
486   // Return: Pointer to the SPD digits array
487
488   // Create an array for 5 chips with 8192 channels each
489   TClonesArray *digitsArray = new TClonesArray("AliITSdigitSPD", kNumberOfChannels);
490
491   return digitsArray;
492 }
493
494
495 //__________________________________________________________________________
496 Bool_t AliITSPreprocessorSPD::CreateGeometryObj(void)
497 {
498   // Creation of the geometry object
499   //
500   // This object is used to get the number of SPD half-staves
501   //
502   // Input : (void)
503   // Output: fGeometryObj
504   // Return: kTRUE if fGeometryObj has been created
505
506   Bool_t status = true;
507
508   // Create geometry object
509   // fGeometryObj = new ...
510   // if (!fGeometryObj) status = kFALSE;
511
512   // Get SPD parameters
513   // fNumberOfColumns = fGeometryObject->GetNumberOfColumns();
514   // fNumberOfRows = fGeometryObject->GetNumberOfRows();
515
516   fNumberOfColumns = kNumberOfColumns;
517   fNumberOfRows = kNumberOfRows;
518
519   return status;
520 }
521
522
523 //__________________________________________________________________________
524 void AliITSPreprocessorSPD::CreateHistograms(void)
525 {
526   // Create the noisy channel histograms
527   //
528   // Input : (void)
529   // Output: Noisy channel histograms
530   // Return: (void)
531
532   TString moduleName;
533   char n[4]; // For storing the module number (maximum module number is 240)
534
535   // Create noisy channel histograms
536   fDigitsHistogram = new TObjArray(fNumberOfModules);
537
538   for (UInt_t i = 0; i < fNumberOfModules; i++)
539     {
540       // Histogram name
541       moduleName = "module_";
542       sprintf(n,"%d",i);
543       moduleName.Append(n);
544
545       (*fDigitsHistogram)[i] = new TH2F(moduleName,"Digits",
546                                         fNumberOfColumns,-.5,(1.*fNumberOfColumns - .5),
547                                         fNumberOfRows,-.5,(1.*fNumberOfRows - .5));
548     }
549 }
550
551
552 //__________________________________________________________________________
553 Bool_t AliITSPreprocessorSPD::FillHistograms(void)
554 {
555   // Fill the histograms with digits (hit positions of unclustered hits)
556   //
557   // (There is one digit histogram per SPD module, i.e. a half-stave, 10 chips)
558   //
559   // Input : No arguments (Empty digit histograms)
560   // Output: Filled digit histograms
561   // Return: kTRUE if histograms are filled with digits, kFALSE otherwise
562
563   AliInfo("Filling noisy channel histograms");
564
565   Bool_t status = kTRUE;
566   AliITSdigitSPD *digitSPD = 0;
567   UInt_t row;
568   UInt_t column;
569   TBranch *digitsSPDBranch;
570   TTree *digitsTree;
571
572   // Get the digits
573   fITSLoader->LoadDigits("read");
574
575   // Create noisy channel histograms
576   AliITSPreprocessorSPD::CreateHistograms();
577
578   // Create an empty SPD digits array
579   TClonesArray *digitsArraySPD = AliITSPreprocessorSPD::CreateDigitsArray();
580
581   // Get number of events
582   UInt_t numberOfEvents = (fRunLoader->TreeE()) ? static_cast<UInt_t>(fRunLoader->TreeE()->GetEntries()) : 0;
583
584   // Make sure we don't try to analyze more data than there actually is
585   if (numberOfEvents < fMaximumNumberOfEvents)
586     {
587       fMaximumNumberOfEvents = numberOfEvents;
588     }
589
590   // Loop over all digits and put them in the corresponding histograms
591   if (numberOfEvents > 0)
592     {
593       for (UInt_t event = 0; event < fMaximumNumberOfEvents; event++)
594         {
595           // Get the current event
596           fRunLoader->GetEvent(event);
597
598           // Get the ITS digits tree
599           digitsTree = fITSLoader->TreeD();
600
601           // Disable all branches except the SPD branch to speed up the reading process
602           digitsTree->SetBranchStatus("ITSDigitsSPD",1);
603           digitsTree->SetBranchStatus("ITSDigitsSDD",0);
604           digitsTree->SetBranchStatus("ITSDigitsSSD",0);
605
606           // Reset the SPD digits array to make sure it is empty
607           digitsArraySPD->Clear();
608
609           // Get the SPD digits branch and set the address
610           digitsSPDBranch = digitsTree->GetBranch("ITSDigitsSPD");
611
612           digitsSPDBranch->SetAddress(&digitsArraySPD);
613
614           if (event%100 == 0) AliInfo(Form("Event #%d", event));
615
616           // Loop over all modules
617           UInt_t numberOfDigits = 0;
618           for (UInt_t module = 0; module < fNumberOfModules; module++)
619             {
620               // Get event data for current module
621               digitsTree->GetEvent(module);
622
623               // Get the number of entries
624               numberOfDigits = digitsArraySPD->GetEntries();
625
626               // Loop over all digits
627               for (UInt_t digit = 0; digit < numberOfDigits; digit++)
628                 {
629                   // Get the current digit
630                   digitSPD = (AliITSdigitSPD*) digitsArraySPD->At(digit);
631                   row = digitSPD->GetCoord1();
632                   column = digitSPD->GetCoord2();
633
634                   // Fill the digits histogram
635                   ((TH2F*)(*fDigitsHistogram)[module])->Fill(column, row);
636
637                 } // end digit loop
638             } // end module loop
639         } // end event loop
640     } // end numberOfEvents > 0
641   else
642     {
643       status = kFALSE;
644     }
645
646   // Cleanup
647   delete digitsArraySPD;
648   digitsArraySPD = 0;
649
650   return status;
651 }
652
653
654 //__________________________________________________________________________
655 Bool_t AliITSPreprocessorSPD::FindDeadChannels(void)
656 {
657   // Locate dead channels
658   //
659   // Input : Filled hit histograms
660   // Output: TObjArray (fBadChannelsObjArray) with all identified dead channels
661   // Return: kTRUE if dead channels have been found
662
663   Bool_t status = kFALSE;
664
665   // Proceed only if properly initialized
666   if (fInit)
667     {
668       if (fVMEMode)
669         {
670           // Initialize counters
671           fNumberOfBadChannels = new Int_t[fNumberOfModules];
672           for (UInt_t module = 0; module < fNumberOfModules; module++)
673             {
674               fNumberOfBadChannels[module] = 0;
675             }
676
677           // Loop over all modules (intentional modularization - DCS will occationally want to
678           // look for noisy channels in certain modules, but not all)
679           fIndex = 0; // Global bad channels array counter (must be reset here)
680           for (UInt_t module = 0; module < fNumberOfModules; module++)
681             {
682               status |= AliITSPreprocessorSPD::FindDeadChannelsInModule(module);
683             }
684         }
685       else
686         {
687           AliError("Dead channels can only be found in data taken with stand-alone VME system");
688         }
689     }
690   else
691     {
692       AliError("Not properly initialized");
693     }
694
695   return status;
696 }
697
698
699 //__________________________________________________________________________
700 Bool_t AliITSPreprocessorSPD::FindDeadChannelsInModule(UInt_t module)
701 {
702   // Locate dead channels
703   // This method assumes that the preprocessor is operator in VME mode.
704   // The algorithm is very simple. It assumes that the data was taken in
705   // a mode where all working SPD pixels should respond as being hit.
706   // fThreshold is used as the limit where everything below this value will
707   // be considered as "dead".
708   //
709   // Input : Filled hit histograms
710   // Output: TObjArray (fBadChannelsObjArray) with all identified noisy channels
711   // Return: kTRUE if dead channels have been found
712
713   // Store the index number for this module
714   fBadChannelsIndexArray[module] = fIndex++;
715
716   UInt_t xBin, numberOfXBins;
717   UInt_t yBin, numberOfYBins;
718   Float_t binContent;
719
720   numberOfXBins = ((TH2F*)(*fDigitsHistogram)[module])->GetNbinsX();
721   numberOfYBins = ((TH2F*)(*fDigitsHistogram)[module])->GetNbinsY();
722
723   // Loop over all bins in this histogram
724   for (xBin = 1; xBin <= numberOfXBins; xBin++)
725     for (yBin = 1; yBin <= numberOfYBins; yBin++)
726       {
727         binContent = ((TH2F*)(*fDigitsHistogram)[module])->GetBinContent(xBin, yBin);
728
729         // Store this channel/bin if outside accepted limits
730         // A channel has to fire MORE times than the fThreshold value, or it will
731         // be considered as "dead"
732         if (binContent <= fThreshold)
733           {
734             // Store the dead channel in the array
735             // The channel object will be deleted in the destructor using the TObjArray Delete() method
736             // (The array will assume ownership of the channel object)
737             AliITSChannelSPD *channel = new AliITSChannelSPD((Int_t)(xBin - 1), (Int_t)(yBin - 1));
738
739             // Store the noisy channel in the array
740             fBadChannelsObjArray->Add(channel);
741
742             // Keep track of the number of bad channels in this module
743             fNumberOfBadChannels[module]++;
744             fIndex += 2;
745
746             // Keep track of the highest module number
747             if (module > fHighestModuleNumber) fHighestModuleNumber = module;
748
749             AliInfo(Form("New dead pixel in (m,c,r) = (%d,%d,%d)", module, xBin - 1, yBin - 1));
750           }
751       } // end bin loop
752
753   return (fNumberOfBadChannels[module] > 0);
754 }
755
756
757 //__________________________________________________________________________
758 Bool_t AliITSPreprocessorSPD::FindNoisyChannels(void)
759 {
760   // Locate noisy channels by searching through the digit histograms
761   // (There is one digit histogram per SPD module, i.e. a half-stave, 5 chips)
762   //
763   // Input : Digits
764   // Output: TObjArray (fBadChannelsObjArray) with all identified noisy channels
765   // Return: kTRUE if noisy channels have been found
766
767   Bool_t status = kFALSE;
768
769   // Proceed only if properly initialized
770   if (fInit)
771     {
772       // (For testing purposes, noisy channels can be inserted here)
773       // SetFakeNoisyChannel(4,10,10);
774
775       // Initialize counters
776       fNumberOfBadChannels = new Int_t[fNumberOfModules];
777       for (UInt_t module = 0; module < fNumberOfModules; module++)
778         {
779           fNumberOfBadChannels[module] = 0;
780         }
781
782       // Scan through all the histogram bins and search for average filling deviations
783       if ((fSelectedAlgorithm == kOptimizedForRealData) || (fSelectedAlgorithm == kOptimizedForRealDataRMS))
784         {
785           // Noisy channel algorithm optimized for real data..........................
786           // Histograms can have any shape (both thresholds and quotients are used)
787           // This algorithm can be used to find noisy channels even if the data was
788           // taken with beam. All channels outside accepted limits (set by fThreshold
789           // and fThresholdRatio) will be marked as noisy
790
791           if (fSelectedAlgorithm == kOptimizedForRealData)
792             {
793               AliInfo("Searching for noisy channels (optimized for real data)");
794             }
795           else
796             {
797               AliInfo("Searching for noisy channels (optimized for real data, RMS version)");
798             }
799
800           // Loop over all modules (intentional modularization - DCS will occationally want to
801           // look for noisy channels in certain modules, but not all)
802           fIndex = 0; // Global bad channels array counter (must be reset here)
803           for (UInt_t module = 0; module < fNumberOfModules; module++)
804             {
805               status |= AliITSPreprocessorSPD::FindNoisyChannelsInModuleAlgo0(module);
806             }
807         } // end algorithm 0
808       else
809         {
810           // Noisy channel algorithm optimized for calibration data...........................
811           // Histograms will/should only contain noisy channels (only thresholds are used)
812           // Calibration histograms should have no background. The calibration data
813           // should have been taken without beam. All channels outside accepted limits
814           // (set by fThreshold) will be marked as noisy
815
816           AliInfo("Searching for noisy channels (optimized for calibration data)");
817
818           // Loop over all modules (intentional modularization - DCS will occationally want to
819           // look for noisy channels in certain modules, but not all)
820           fIndex = 0; // Global bad channels array counter (must be reset here)
821           for (UInt_t module = 0; module < fNumberOfModules; module++)
822             {
823               status |= AliITSPreprocessorSPD::FindNoisyChannelsInModuleAlgo1(module);
824             }
825         } // end algorithm 1
826     } // end if fInit
827   else
828     {
829       AliError("Not properly initialized");
830     }
831
832   return status;
833 }
834
835
836 //__________________________________________________________________________
837 Bool_t AliITSPreprocessorSPD::FindNoisyChannelsInModuleAlgo0(UInt_t module)
838 {
839   // Locate the noisy channels in a module (optimized for real data)
840   //
841   // For each channel in these histograms, the algorithm checks the average
842   // filling of the neighboring 3, 5 or 8 channels (depending on the location
843   // of the channel in question; corner, border or inside), or compares with the
844   // RMS of the neighbors. If the average is "much" lower, the channel will be
845   // considered to be noisy. The default noisy-to-normal fraction is stored in the
846   // fThresholdRatio varible. It can be set with the SetThresholdRatio method.
847   // The channel also has to have fired more than a certain minimum, fThreshold.
848   // It can be set with the SetThreshold method.
849   //
850   // To avoid difficulties with noisy channels that occur in pairs, the
851   // neighboring channel with largest number of fillings will be removed from
852   // the average calculation.
853   //
854   // NOTE: Since this method modifies the fBadChannelsObjArray and fBadChannelsIndexArray
855   //       it is essential to initialize the fIndex counter before calling this module
856   //       the first time. The bad channel data does not have to be ordered per module
857   //       in the fBadChannelsObjArray, but the indices of where the data of a certain module
858   //       starts has to be correct. A wrong fIndex can lead to segmentation violation
859   //
860   // Input : module number, filled digit histograms
861   // Output: TObjArray (fBadChannelsObjArray) with all identified noisy channels,
862   //         Int_t[] (fBadChannelsIndexArray) with fBadChannelsObjArray module indices,
863   //         number of noisy channels in this module (global variable fNumberOfBadChannels[module])
864   // Return: kTRUE if there are noisy channels in this module
865
866   // Store the index number for this module
867   fBadChannelsIndexArray[module] = fIndex++;
868
869   UInt_t xBin, numberOfXBins;
870   UInt_t yBin, numberOfYBins;
871   UInt_t neighborXBin;
872   UInt_t neighborYBin;
873   UInt_t numberOfNeighboringBins;
874   Float_t binContent;
875   Float_t sumBinContent;
876   Float_t neighborBinContent;
877   Float_t maxBinContent;
878   Float_t averageBinContent;
879   Float_t ratio;
880
881   numberOfXBins = (UInt_t)((TH2F*)(*fDigitsHistogram)[module])->GetNbinsX();
882   numberOfYBins = (UInt_t)((TH2F*)(*fDigitsHistogram)[module])->GetNbinsY();
883
884   // Loop over all bins in this histogram
885   for (xBin = 1; xBin <= numberOfXBins; xBin++)
886     for (yBin = 1; yBin <= numberOfYBins; yBin++)
887       {
888         numberOfNeighboringBins = 0;
889         averageBinContent = 0.;
890         sumBinContent = 0.;
891         binContent = ((TH2F*)(*fDigitsHistogram)[module])->GetBinContent(xBin, yBin);
892         maxBinContent = 0.;
893
894         // Calculate the average pixel level on the surrounding pixels
895         for (neighborXBin = xBin - 1; neighborXBin <= xBin + 1; neighborXBin++)
896           for (neighborYBin = yBin - 1; neighborYBin <= yBin + 1; neighborYBin++)
897             {
898               if ( !((neighborXBin == xBin) && (neighborYBin == yBin)) )
899                 {
900                   // Only update the number of neighboring bins when we are not on the border
901                   if ((neighborXBin > 0) && (neighborXBin <= numberOfXBins+1) &&
902                       (neighborYBin > 0) && (neighborYBin <= numberOfYBins+1))
903                     {
904                       neighborBinContent =
905                         ((TH2F*)(*fDigitsHistogram)[module])->GetBinContent(neighborXBin, neighborYBin);
906
907                       if (fSelectedAlgorithm == kOptimizedForRealDataRMS)
908                         {
909                           // RMS
910                           sumBinContent += neighborBinContent*neighborBinContent;
911                         }
912                       else
913                         {
914                           // Geometrical mean
915                           sumBinContent += neighborBinContent;
916                         }
917
918                       if (neighborBinContent > maxBinContent) maxBinContent = neighborBinContent;
919
920                       numberOfNeighboringBins++;
921                     }
922                 }
923             }
924
925         // Calculate the average. Remove the largest neighboring bin
926         // (Correction for potential clusters of noisy channels)
927         if (fSelectedAlgorithm == kOptimizedForRealDataRMS)
928           {
929             // Square the max bin content before removing it from the average calculation
930             maxBinContent *= maxBinContent;
931
932             // RMS
933             averageBinContent = sqrt((sumBinContent - maxBinContent)/(Float_t)(numberOfNeighboringBins - 1));
934           }
935         else
936           {
937             // Geometrical mean
938             averageBinContent = (sumBinContent - maxBinContent)/(Float_t)(numberOfNeighboringBins - 1);
939           }
940
941         // Store this channel/bin if outside accepted limits
942         // The threshold ratio is the threshold for the current bin content divided by the
943         // average neighboring bin contents. The threshold bin content is the minimum number of
944         // times a channel has to have fired to be called noisy
945         ratio = (averageBinContent > 0) ? binContent/averageBinContent : 0.;
946         if ( ((ratio >= fThresholdRatio) || (ratio == 0.)) && (binContent >= fThreshold) )
947           {
948             // Store the noisy channel in the array
949             // The channel object will be deleted in the destructor using the TObjArray Delete() method
950             // (The array will assume ownership of the channel object)
951             AliITSChannelSPD *channel = new AliITSChannelSPD((Int_t)(xBin - 1), (Int_t)(yBin - 1));
952
953             // Store the noisy channel in the array
954             fBadChannelsObjArray->Add(channel);
955
956             // Keep track of the number of bad channels in this module
957             fNumberOfBadChannels[module]++;
958             fIndex += 2;
959
960             // Keep track of the highest module number
961             if (module > fHighestModuleNumber) fHighestModuleNumber = module;
962
963             AliInfo(Form("New noisy pixel in (m,c,r) = (%d,%d,%d)", module, xBin - 1, yBin - 1));
964             AliInfo(Form("- Noisy pixel fired %d times, average neighborhood: %f",(Int_t)binContent,averageBinContent));
965           }
966       } // end bin loop
967
968   return (fNumberOfBadChannels[module] > 0);
969 }
970
971
972 //__________________________________________________________________________
973 Bool_t AliITSPreprocessorSPD::FindNoisyChannelsInModuleAlgo1(UInt_t module)
974 {
975   // Locate the noisy channels in a module (optimized for calibration data)
976   //
977   // This algorithm locates noisy channels by assuming original data was taken
978   // in calibration mode. This should be done without beam and will thus only
979   // contain data corresponding to background and noisy channels. The latter
980   // should be clearly visible in this data so this algorithm simply assumes
981   // that all histogram bins that are filled more than fThreshold times are
982   // noisy.
983   //
984   // NOTE: Since this method modifies the fBadChannelsObjArray and fBadChannelsIndexArray
985   //       it is essential to initialize the fIndex counter before calling this module
986   //       the first time. The bad channel data does not have to be ordered per module
987   //       in the fBadChannelsObjArray, but the indices of where the data of a certain module
988   //       starts has to be correct. A wrong fIndex can lead to segmentation violation
989   //
990   // Input : module number, filled digit histograms
991   // Output: TObjArray (fBadChannelsObjArray) with all identified noisy channels,
992   //         Int_t[] (fBadChannelsIndexArray) with fBadChannelsObjArray module indices,
993   //         number of noisy channels in this module (global variable fNumberOfBadChannels[module])
994   // Return: kTRUE if there are noisy channels in this module
995
996   // Store the index number for this module
997   fBadChannelsIndexArray[module] = fIndex++;
998
999   UInt_t xBin, numberOfXBins;
1000   UInt_t yBin, numberOfYBins;
1001   Float_t binContent;
1002
1003   numberOfXBins = ((TH2F*)(*fDigitsHistogram)[module])->GetNbinsX();
1004   numberOfYBins = ((TH2F*)(*fDigitsHistogram)[module])->GetNbinsY();
1005
1006   // Loop over all bins in this histogram
1007   for (xBin = 1; xBin <= numberOfXBins; xBin++)
1008     for (yBin = 1; yBin <= numberOfYBins; yBin++)
1009       {
1010         binContent = ((TH2F*)(*fDigitsHistogram)[module])->GetBinContent(xBin, yBin);
1011
1012         // Store this channel/bin if outside accepted limits
1013         // The threshold bin content is the minimum number of times a channel has to have
1014         // fired to be called noisy
1015         if (binContent >= fThreshold)
1016           {
1017             // Store the noisy channel in the array
1018             // The channel object will be deleted in the destructor using the TObjArray Delete() method
1019             // (The array will assume ownership of the channel object)
1020             AliITSChannelSPD *channel = new AliITSChannelSPD((Int_t)(xBin - 1), (Int_t)(yBin - 1));
1021
1022             // Store the noisy channel in the array
1023             fBadChannelsObjArray->Add(channel);
1024
1025             // Keep track of the number of bad channels in this module
1026             fNumberOfBadChannels[module]++;
1027             fIndex += 2;
1028
1029             // Keep track of the highest module number
1030             if (module > fHighestModuleNumber) fHighestModuleNumber = module;
1031
1032             AliInfo(Form("New noisy pixel in (m,c,r) = (%d,%d,%d)", module, xBin - 1, yBin - 1));
1033             AliInfo(Form("- Noisy pixel fired %d times",(Int_t)binContent));
1034           }
1035       } // end bin loop
1036
1037   return (fNumberOfBadChannels[module] > 0);
1038 }
1039
1040
1041 //__________________________________________________________________________
1042 void AliITSPreprocessorSPD::PrintChannels(void)
1043 {
1044   // Print all found bad channels to stdout
1045   //
1046   // Input : fBadChannelsObjArray
1047   // Output: (dump to stdout)
1048   // Return: (void)
1049
1050   Int_t i = 0;
1051   Int_t j = 0;
1052   AliITSChannelSPD *channel = 0;
1053
1054   // Print the bad channels stores in the array
1055   AliInfo("\nModule #\tColumn #\tRow #\n------------------------------------------------");
1056   for (UInt_t module = 0; module < fNumberOfModules; module++)
1057     {
1058       j = 0;
1059       while (j < fNumberOfBadChannels[module])
1060         {
1061           channel = (AliITSChannelSPD *) fBadChannelsObjArray->At(i++);
1062           std::cout << module << "\t\t" << channel->GetColumn() << "\t\t" << channel->GetRow() << std::endl;
1063
1064           // Go to next bad channel
1065           j++;
1066         }
1067     }
1068
1069   AliInfo(Form("%d bad channels were found", fBadChannelsObjArray->GetEntries()));
1070 }
1071
1072
1073 //__________________________________________________________________________
1074 void AliITSPreprocessorSPD::MarkNoisyChannels(void)
1075 {
1076   // WARNING: THIS METHOD DOESN'T WORK!!!
1077   //
1078   // Mark all identified noisy channels
1079   //
1080   // Input : List of noisy channels, original digits tree
1081   // Output: New digits tree containing SPD digits marked when noisy
1082   // Return: (void)
1083   //
1084   // The original digits tree (digitsTree) is cloned except for the SPD branch (ITSDigitSPD).
1085   // This branch is then redefined for each event and will contain all the original
1086   // information. All known noisy channels will be marked by using the TObject status bits
1087   // according to the following scheme. Dead channels are included for completeness. Note
1088   // that a dead channel will NEVER show up among digits..
1089   //
1090   // Interpretation of digit status bits (LSB):
1091   // Dead channel   Noisy channel   |  Integer
1092   // -----------------------------------------
1093   //      0               0         |     0
1094   //      0               1         |     1
1095   //      1               0         |     2
1096   //
1097   // meaning e.g. that a channel that is noisy will have the first bit set in its status bits
1098
1099   // Do not continue unless we are processing DAQ data
1100   if (!fVMEMode)
1101     {
1102       AliInfo("Marking bad channels");
1103
1104       // Create the storage container that will be used to access the bad channels
1105       if (!fBadChannelsContainer)
1106         {
1107           // Add the bad channels array to the storage container
1108           // (ownership is passed to the AliRunDataStorage object)
1109           fBadChannelsContainer = new AliITSBadChannelsSPD();
1110
1111           // Convert the bad channels from TObjArray to Int_t[]
1112           AliITSPreprocessorSPD::ConvertObjToIntArray();
1113
1114           // Store the arrays in the bad channels container object
1115           const Int_t kBadChannelsArraySize =
1116             2*fBadChannelsObjArray->GetEntries() + fNumberOfModules;
1117           fBadChannelsContainer->Put(fBadChannelsIntArray, kBadChannelsArraySize,
1118                                      fBadChannelsIndexArray, fNumberOfModules);
1119         }
1120
1121       // Create the bad channels helper object
1122       // (will be used to find a bad channel within a TObjArray)
1123       AliITSBadChannelsAuxSPD *aux = new AliITSBadChannelsAuxSPD();
1124
1125       AliITSdigitSPD *digitSPD = 0;
1126       UInt_t numberOfDigits;
1127       Int_t newDigit[3];
1128       Bool_t mark = kFALSE;
1129
1130       TBranch *digitsBranch = 0;
1131       TTree *digitsTree;
1132
1133       // Create an empty SPD digit array
1134       TObjArray *digitsArraySPD = new TObjArray();
1135
1136       // Get the digits in update mode (we want to modify them if there are noisy channels)
1137       fITSLoader->UnloadDigits();
1138       fITSLoader->LoadDigits("update");
1139
1140       // Get the number of events
1141       UInt_t numberOfEvents = (fRunLoader->TreeE()) ? static_cast<UInt_t>(fRunLoader->TreeE()->GetEntries()) : 0;
1142
1143       // Loop over all events
1144       for (UInt_t event = 0; event < numberOfEvents; event++)
1145         {
1146           if (event%100 == 0) AliInfo(Form("Event #%d", event));
1147
1148           // Get the current event
1149           fRunLoader->GetEvent(event);
1150
1151           // Get the ITS digits tree
1152           digitsTree = fITSLoader->TreeD();
1153
1154           // Get SPD branch that will contain all digits with marked noisy channels
1155           digitsBranch = digitsTree->GetBranch("ITSDigitsSPD");
1156           digitsBranch->SetAddress(&digitsArraySPD);
1157
1158           // Get the stored number of modules
1159           UInt_t numberOfModules = (Int_t)digitsTree->GetEntries();
1160           TObjArray **newDigitsArraySPD = new TObjArray*[numberOfModules];
1161
1162           Int_t *digitNumber = new Int_t[numberOfModules];
1163           for (UInt_t m = 0; m < numberOfModules; m++)
1164             {
1165               newDigitsArraySPD[m] = new TObjArray();
1166               digitNumber[m] = 0;
1167             }
1168
1169           AliInfo(Form("ent = %d", (Int_t)digitsTree->GetEntries()));
1170
1171           // Reset the SPD digit arrays to make sure they are empty
1172           digitsArraySPD->Clear();
1173
1174           // Get the SPD digits branch from the original digits tree and set the address
1175           digitsBranch = digitsTree->GetBranch("ITSDigitsSPD");
1176           digitsBranch->SetAddress(&digitsArraySPD);
1177
1178           // Loop over all modules
1179           for (UInt_t module = 0; module < fNumberOfModules; module++)
1180             {
1181               // Get event data for current module
1182               digitsTree->GetEvent(module);
1183
1184               // Get the hits in the current module
1185               TObjArray *moduleObjArray = fBadChannelsContainer->CreateModuleObjArray(module);
1186
1187               // Get the number of entries
1188               numberOfDigits = digitsArraySPD->GetEntries();
1189
1190               // Loop over all digits and all channels
1191               for (UInt_t digit = 0; digit < numberOfDigits; digit++)
1192                 {
1193                   // Get the current digit
1194                   digitSPD = (AliITSdigitSPD*) digitsArraySPD->At(digit);
1195                   newDigit[0] = digitSPD->GetCoord1(); // row
1196                   newDigit[1] = digitSPD->GetCoord2(); // column
1197                   newDigit[2] = digitSPD->GetSignal(); // signal
1198
1199                   // Check if this channel is noisy
1200                   // (Compare with all stored channels in the bad channels array)
1201                   if (aux->Find(digitSPD, moduleObjArray))
1202                     {
1203                       // Set the mark flag and break the loop
1204                       mark = kTRUE;
1205                     }
1206
1207                   // Store this digit in the SPD digits array using a placement new operation
1208                   new ((*newDigitsArraySPD[module])[digitNumber[module]]) AliITSdigitSPD(newDigit);
1209
1210                   // Mark it if noisy and store in the noisy channel array
1211                   if (mark)
1212                     {
1213                       // Store this digit in the marked SPD digits array using a placement new operation
1214                       //new ((*badChannels[m])[numberOfBadChannels[m]]) AliITSChannelSPD(newBadChannel);
1215                       //new ((*newDigitsArraySPD[module])[digitNumber[module]]) AliITSdigitSPD(newDigit);
1216
1217                       // Mark the original channel as noisy
1218                       ((*newDigitsArraySPD[module])[digitNumber[module]])->SetBit(kNoisyChannel);
1219
1220                       mark = kFALSE;
1221                     }
1222
1223                   digitNumber[module]++;
1224
1225                 } // end digit loop
1226
1227               // Cleanup
1228               delete moduleObjArray;
1229               moduleObjArray = 0;
1230
1231             } // end module loop
1232
1233           digitsBranch->Reset();
1234           digitsBranch->ResetAddress();
1235
1236           // Cleanup
1237           delete digitsArraySPD;
1238           digitsArraySPD = 0;
1239           digitsTree->Reset();
1240
1241           // WHY THIS RANGE?????????????????????????????????????????????????????????????????????
1242           for (UInt_t n = 0; n < event; n++)
1243             {
1244               digitsTree->SetBranchAddress("ITSDigitsSPD", &newDigitsArraySPD[n]);
1245               digitsTree->Fill();
1246             }
1247
1248           digitsTree->AutoSave();
1249
1250           // Cleanup
1251           for (UInt_t n = 0; n < event; n++)
1252             {
1253               delete newDigitsArraySPD[n];
1254             }
1255           delete [] newDigitsArraySPD;
1256           newDigitsArraySPD = 0;
1257           delete [] digitNumber;
1258           digitNumber = 0;
1259           delete digitsTree;
1260           digitsTree = 0;
1261
1262         } // end loop over all events
1263
1264       // Unload the digits
1265       fITSLoader->UnloadDigits();
1266
1267       // Cleanup
1268       delete aux;
1269       aux = 0;
1270     }
1271 }
1272
1273
1274
1275 //__________________________________________________________________________
1276 Bool_t AliITSPreprocessorSPD::Store(AliCDBId &id, AliCDBMetaData *md)
1277 {
1278   // Store the bad channels object in the calibration database
1279   // (See the corresponding run macro for further explanations)
1280   //
1281   // Input : fBadChannelsObjArray (now containing all found bad channels), object meta data
1282   // Output: Database file containing the bad channels
1283   // Return: kTRUE if successful
1284
1285   Bool_t status = kFALSE;
1286
1287   AliInfo("Storing bad channels");
1288
1289   // Add the bad channels array to the storage container
1290   // (ownership is passed to the AliRunDataStorage object)
1291   fBadChannelsContainer = new AliITSBadChannelsSPD();
1292
1293   // Convert the bad channels from TObjArray to Int_t[]
1294   AliITSPreprocessorSPD::ConvertObjToIntArray();
1295
1296   // Store the arrays in the bad channels container object
1297   const Int_t kBadChannelsArraySize =
1298     2*fBadChannelsObjArray->GetEntries() + fNumberOfModules;
1299   fBadChannelsContainer->Put(fBadChannelsIntArray, kBadChannelsArraySize,
1300                              fBadChannelsIndexArray, fNumberOfModules);
1301
1302   // Store the container
1303  if(!AliCDBManager::Instance()->IsDefaultStorageSet()) {
1304         AliError("No storage set!");
1305         return status;
1306   }
1307   
1308   if (AliCDBManager::Instance()->GetDefaultStorage()->Put(fBadChannelsContainer, id, md))
1309     {
1310       AliInfo("Bad channels object stored in database");
1311       status = kTRUE;
1312     }
1313   else
1314     {
1315       AliError("Failed to store object in database");
1316     }
1317
1318   return status;
1319 }
1320
1321
1322 //__________________________________________________________________________
1323 void AliITSPreprocessorSPD::ConvertObjToIntArray(void)
1324 {
1325   // Convert the bad channel TObjArray to an Int_t array
1326   //
1327   // Input : fBadChannelsObjArray (now containing all found bad channels)
1328   // Output: fBadChannelsIntArray
1329   // Return: (void)
1330   //
1331   // Data encoding:
1332   //                The TObjArray of this class (fBadChannelsObjArray) is converted to a sequential
1333   //                Int_t array (fBadChannelsIntArray) in this method. For each module, the first
1334   //                stored number is the number of bad channels in the current module. This is
1335   //                followed by all the columns and rows of the bad channels:
1336   //
1337   //                badChannelsArray =
1338   //                | N(m) | col0 | row0 | .. | colN(m) | N(m+1) | col0 | row0 | ...
1339   //                .   .......... module m .........   .   .... module m+1 ......
1340   //
1341   //                The bad channels index array (fBadChannelsIndexArray) contains the indices of
1342   //                the badChannelsArray, i.e. where the bad channels in certain module starts:
1343   //
1344   //                fBadChannelsObjArray =
1345   //                | i0 | i1 | .. | iM | (where M = the number of SPD modules)
1346   //
1347   //                e.g. i1 corresponds to the index of the badChannelsArray where N(1) is stored,
1348   //                i.e. the number of bad channels for module 1
1349
1350   const Int_t kBadChannelsArraySize =
1351     2*fBadChannelsObjArray->GetEntries() + fNumberOfModules;
1352   fBadChannelsIntArray = new Int_t[kBadChannelsArraySize]; // Will be deleted in dtor
1353   AliITSChannelSPD *channel = 0;
1354   Int_t i = 0;
1355   Int_t j = 0;
1356   Int_t k = 0;
1357
1358   // Loop over all modules
1359   for (UInt_t module = 0; module < fNumberOfModules; module++)
1360     {
1361       // Encode the number of bad channels of the current module
1362       fBadChannelsIntArray[k++] = fNumberOfBadChannels[module];
1363
1364       // The columns and rows of the fBadChannelsObjArray will be stored sequentially
1365       // in the Int_t array
1366       j = 0;
1367       while (j < fNumberOfBadChannels[module])
1368         {
1369           channel = (AliITSChannelSPD *) fBadChannelsObjArray->At(i++);
1370           fBadChannelsIntArray[k++] = channel->GetColumn();
1371           fBadChannelsIntArray[k++] = channel->GetRow();
1372
1373           // Go to next bad channel
1374           j++;
1375         }
1376     }
1377 }