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