1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.2 2006/02/03 11:31:18 masera
19 Calibration framework improved (E. Crescio)
21 Revision 1.1 2005/10/11 12:31:50 masera
22 Preprocessor classes for SPD (Paul Nilsson)
26 ///////////////////////////////////////////////////////////////////////////
27 // AliITSPreprocessorSPD implementation by P. Nilsson 2005
28 // MAIN AUTHOR/CONTACT: Paul.Nilsson@cern.ch
30 // The purpose of this class is to provide algorithms for identification
31 // of "bad channels" such as dead channels and noisy channels in the SPD
33 // Examples on how to use this class can be found in the macros:
35 // .findNoisyChannels.C (Locate and store noisy channels in the CDB)
36 // .readNoisyChannels.C (Read noisy channels from the CDB)
37 // .findDeadChannels.C (Locate and store dead channels in the CDB)
38 // .readDeadChannels.C (Read dead channels from the CDB)
40 // Modified by D. Elia, H. Tydesjo
41 // March 2006: Mixed up coordinates, bug fixed
43 ///////////////////////////////////////////////////////////////////////////
46 #include "AliITSPreprocessorSPD.h"
47 #include "AliCDBEntry.h"
48 #include "AliITSCalibrationSPD.h"
49 ClassImp(AliITSPreprocessorSPD)
52 //__________________________________________________________________________
53 AliITSPreprocessorSPD::AliITSPreprocessorSPD(void):
58 fMaximumNumberOfEvents(1000000),
59 fHighestModuleNumber(0),
60 fSelectedAlgorithm(kOptimizedForRealData),
61 fGeometryMode(kALICEGeometry),
62 fNumberOfBadChannels(0),
66 fBadChannelsObjArray(0),
67 fBadChannelsIntArray(0),
68 fBadChannelsIndexArray(0),
69 fBadChannelsContainer(0)
71 // Default constructor for the SPD preprocessor
73 // Initialization has to be done by hand using Set* methods and the Open method
77 // Return: <empty/uninitialized AliITSPreprocessorSPD object>
81 //__________________________________________________________________________
82 AliITSPreprocessorSPD::AliITSPreprocessorSPD(const char *fileName, const char *mode,
83 const char *fileNameg, const Int_t maxNumberOfEvents):
88 // Standard constructor for the SPD preprocessor
90 // Input : Name of digit file
92 // Return: Initialized SPD preprocessor object
95 AliITSPreprocessorSPD::Init();
96 AliITSPreprocessorSPD::SetMaximumNumberOfEvents(maxNumberOfEvents);
98 // Open and read the galice and digit files
99 if (!AliITSPreprocessorSPD::Open(fileName, mode, fileNameg))
101 AliError("Failed to open file");
107 //__________________________________________________________________________
108 AliITSPreprocessorSPD::AliITSPreprocessorSPD(const AliITSPreprocessorSPD &prep) :
111 // Default copy constructor
112 // Notice that only pointer addresses are copied!
113 // Memory allocations of new objects are not done.
115 fITSLoader = prep.fITSLoader;
116 fRunLoader = prep.fRunLoader;
117 fThresholdRatio = prep.fThresholdRatio;
118 fThreshold = prep.fThreshold;
119 fMaximumNumberOfEvents = prep.fMaximumNumberOfEvents;
120 fHighestModuleNumber = prep.fHighestModuleNumber;
121 fSelectedAlgorithm = prep.fSelectedAlgorithm;
122 fGeometryMode = prep.fGeometryMode;
123 fNumberOfBadChannels = prep.fNumberOfBadChannels;
125 fVMEMode = prep.fVMEMode;
126 fDigitsHistogram = prep.fDigitsHistogram;
127 fBadChannelsObjArray = prep.fBadChannelsObjArray;
128 fBadChannelsIntArray = prep.fBadChannelsIntArray;
129 fBadChannelsIndexArray = prep.fBadChannelsIndexArray;
130 fBadChannelsContainer = prep.fBadChannelsContainer;
134 //__________________________________________________________________________
135 AliITSPreprocessorSPD& AliITSPreprocessorSPD::operator=(const AliITSPreprocessorSPD &prep)
137 // Assignment operator
139 AliError("Not implemented");
141 if (this != &prep) { } // Do not delete this line
147 //__________________________________________________________________________
148 AliITSPreprocessorSPD::~AliITSPreprocessorSPD(void)
150 // Destructor for the SPD preprocessor
152 if (fDigitsHistogram)
154 // try fDigitsHistogram->Delete(); if the following crashes
155 for (UInt_t module = 0; module < fNumberOfModules; module++)
157 (*fDigitsHistogram)[module]->Delete();
159 delete fDigitsHistogram;
160 fDigitsHistogram = 0;
163 if (fNumberOfBadChannels)
165 delete [] fNumberOfBadChannels;
166 fNumberOfBadChannels = 0;
169 if (fBadChannelsIntArray)
171 delete [] fBadChannelsIntArray;
172 fBadChannelsIntArray = 0;
175 if (fBadChannelsIndexArray)
177 delete [] fBadChannelsIndexArray;
178 fBadChannelsIndexArray = 0;
181 if (fBadChannelsObjArray)
183 fBadChannelsObjArray->Delete();
184 delete fBadChannelsObjArray;
185 fBadChannelsObjArray = 0;
191 delete fBadChannelsContainer;
192 fBadChannelsContainer = 0;
196 //__________________________________________________________________________
197 void AliITSPreprocessorSPD::Init(void)
199 // Initialization of the SPD preprocessor
202 // Output: Several logistical variables
205 // Default maximum number of events per histogram
206 fMaximumNumberOfEvents = 1000000;
208 // Default noisy channel removal algorithm
209 fSelectedAlgorithm = kOptimizedForRealData;
211 // Default noisy channel threshold and threshold ratio
212 // (threshold for current bin content divided by the average neighboring bin contents)
214 fThresholdRatio = 5.;
216 // Set default geometry mode (ALICE geometry). This also sets fNumberOfModules
217 AliITSPreprocessorSPD::SetGeometryMode(kALICEGeometry);
219 // The highest module number with found bad channels
220 fHighestModuleNumber = 0;
222 // Assume input is not from a VME file
225 // Initialization is complete
230 //__________________________________________________________________________
231 void AliITSPreprocessorSPD::SetGeometryMode(UInt_t mode)
233 // Set the geometry mode
235 // Input : Geometry mode (either kTestBeamGeometry or kALICEGeometry)
239 fGeometryMode = mode;
241 // In case of an input VME file, the number of modules has already been fixed.
242 // Do not try to change it
245 if (mode == kALICEGeometry)
247 fNumberOfModules = kNumberOfSPDModules;
249 else if (mode == kTestBeamGeometry)
251 fNumberOfModules = kNumberOfTestBeamSPDModules;
255 AliError("Unknown geometry mode, defaults to ALICE geometry");
256 fNumberOfModules = kNumberOfSPDModules;
262 //__________________________________________________________________________
263 void AliITSPreprocessorSPD::SetFakeNoisyChannel(Int_t module, Int_t column, Int_t row)
265 // Introduce a fake noisy channel in the hit histograms
267 // Input : Module, column and row numbers
268 // Output: Updated hit histograms
271 if ((UInt_t)module < fNumberOfModules)
273 ((TH2F*)(*fDigitsHistogram)[module])->Fill(column, row, 1000000);
277 AliError("No such module number");
282 //__________________________________________________________________________
283 Bool_t AliITSPreprocessorSPD::Open(const char *fileName, const char *mode, const char *fileNameg)
287 // Input : Name and mode of ITS digit file, name of galice file
289 // Return: kTRUE if loaders succeed
291 Bool_t status = kFALSE;
292 Bool_t status0 = kFALSE;
293 Bool_t status1 = kFALSE;
294 Bool_t status2 = kFALSE;
295 Bool_t status3 = kFALSE;
297 // Only proceed if initialization has been done
300 TString m = (TString)mode;
301 if (m == "daq" || m == "DAQ")
303 // Open the data file and get the run loader
304 fRunLoader = AliRunLoader::Open(fileNameg);
307 // Get the gAlice object
308 status0 = AliITSPreprocessorSPD::GetgAlice();
310 // Get the ITS digits
311 if (status0) status1 = AliITSPreprocessorSPD::GetITSDigits(fileName);
313 // Create the test beam object
314 if (status1) status2 = AliITSPreprocessorSPD::CreateGeometryObj();
316 // Fill histograms with DAQ digits
317 if (status2) status3 = AliITSPreprocessorSPD::FillHistograms();
319 status = status0 & status1 & status2 & status3;
323 AliError("Failed to get the run loader");
326 else if (m == "vme" || m == "VME")
328 // Open the VME file that contains histograms with fired channels as read by the stand-alone VME system
329 TFile *vmeFile = TFile::Open(fileName);
331 if (!vmeFile->IsOpen())
333 AliError("Could not open VME input file");
337 // Get the histograms from the VME file that contains all fired channels
338 status0 = AliITSPreprocessorSPD::GetVMEHistograms(vmeFile);
340 // Create the test beam object
341 if (status0) status1 = AliITSPreprocessorSPD::CreateGeometryObj();
344 // This boolean will be used to override any attempts of changing the number of modules by the user
345 // with the SetGeometryMode method. For VME files the number of modules will entirely be determined by
346 // the input VME root file, i.e. by the number of histograms in this file
349 status = status0 & status1;
353 AliError("Unknown filetype - assuming DAQ file");
356 // At this stage, the final number of modules will be known (safe to define arrays)
357 // In case data is read from a VME root file, it will not be known before
360 // Create the arrays for bookkeeping and storing the noisy channels
361 if (!fBadChannelsObjArray)
363 fBadChannelsObjArray = new TObjArray();
365 if (!fBadChannelsIndexArray)
367 // This array will contain the begin and end indices for each module, i.e. where to begin
368 // and stop reading the fBadChannelsObjArray for a certain module.
369 // The last position of the array will contain the end index of the last module
370 fBadChannelsIndexArray = new Int_t[fNumberOfModules + 1];
371 for (UInt_t module = 0; module < fNumberOfModules + 1; module++)
373 fBadChannelsIndexArray[module] = 0;
380 AliError("SPD preprocessor not initialized. Can't load digits");
387 //__________________________________________________________________________
388 Bool_t AliITSPreprocessorSPD::GetgAlice(void)
390 // Get the gAlice object
393 // Output: gAlice object
394 // Return: kTRUE if the gAlice object was found
396 Bool_t status = kTRUE;
399 fRunLoader->LoadgAlice();
400 if (!fRunLoader->GetAliRun())
402 AliError("gAlice not found on file. Aborting.");
410 //__________________________________________________________________________
411 Bool_t AliITSPreprocessorSPD::GetVMEHistograms(TFile *vmeFile)
413 // Get pre-filled digit histograms from input VME file
415 // Input : pointer to VME file
417 // Return: kTRUE if histograms are found on file
419 Bool_t status = kFALSE;
421 // Get the file directory
422 TDirectory *dir = (TDirectory *)vmeFile;
424 // Get the number of keys (histograms in this case corresponding to modules/ladders)
425 fNumberOfModules = dir->GetNkeys();
426 if ((fNumberOfModules > 0) && (fNumberOfModules <= kNumberOfSPDModules))
430 // Create bad channel histograms
431 fDigitsHistogram = new TObjArray(fNumberOfModules);
433 // Create a key iterator
434 TIter nextkey(dir->GetListOfKeys());
437 // Loop over all objects, read them in to memory one by one
439 while ( (key = (TKey *)nextkey()) )
441 (*fDigitsHistogram)[module++] = (TH2F *)key->ReadObj();
444 if (module > fNumberOfModules)
446 AliError("Wrong number of keys in VME file");
454 AliError("Wrong number of histograms in VME file");
459 //__________________________________________________________________________
460 Bool_t AliITSPreprocessorSPD::GetITSDigits(const char *fileName)
462 // Get the ITS digits
464 // Input : name of digits file
465 // Output: fITSLoader object, ITS digits
466 // Return: kTRUE if loader succeed
468 Bool_t status = kTRUE;
470 // Load the gAlice and the header
471 fRunLoader->LoadgAlice();
472 fRunLoader->LoadHeader();
474 // Get the ITS loader
475 fITSLoader = (AliITSLoader*) fRunLoader->GetLoader("ITSLoader");
478 AliError("ITS loader not found");
483 // Open the digits file
484 fITSLoader->SetDigitsFileName(fileName);
491 //__________________________________________________________________________
492 TClonesArray* AliITSPreprocessorSPD::CreateDigitsArray(void) const
494 // Creation of the digits array
498 // Return: Pointer to the SPD digits array
500 // Create an array for 5 chips with 8192 channels each
501 TClonesArray *digitsArray = new TClonesArray("AliITSdigitSPD", kNumberOfChannels);
507 //__________________________________________________________________________
508 Bool_t AliITSPreprocessorSPD::CreateGeometryObj(void)
510 // Creation of the geometry object
512 // This object is used to get the number of SPD half-staves
515 // Output: fGeometryObj
516 // Return: kTRUE if fGeometryObj has been created
518 Bool_t status = true;
520 // Create geometry object
521 // fGeometryObj = new ...
522 // if (!fGeometryObj) status = kFALSE;
524 // Get SPD parameters
525 // fNumberOfColumns = fGeometryObject->GetNumberOfColumns();
526 // fNumberOfRows = fGeometryObject->GetNumberOfRows();
528 fNumberOfColumns = kNumberOfColumns;
529 fNumberOfRows = kNumberOfRows;
535 //__________________________________________________________________________
536 void AliITSPreprocessorSPD::CreateHistograms(void)
538 // Create the noisy channel histograms
541 // Output: Noisy channel histograms
545 char n[4]; // For storing the module number (maximum module number is 240)
547 // Create noisy channel histograms
548 fDigitsHistogram = new TObjArray(fNumberOfModules);
550 for (UInt_t i = 0; i < fNumberOfModules; i++)
553 moduleName = "module_";
555 moduleName.Append(n);
557 (*fDigitsHistogram)[i] = new TH2F(moduleName,"Digits",
558 fNumberOfColumns,-.5,(1.*fNumberOfColumns - .5),
559 fNumberOfRows,-.5,(1.*fNumberOfRows - .5));
564 //__________________________________________________________________________
565 Bool_t AliITSPreprocessorSPD::FillHistograms(void)
567 // Fill the histograms with digits (hit positions of unclustered hits)
569 // (There is one digit histogram per SPD module, i.e. a half-stave, 10 chips)
571 // Input : No arguments (Empty digit histograms)
572 // Output: Filled digit histograms
573 // Return: kTRUE if histograms are filled with digits, kFALSE otherwise
575 AliInfo("Filling noisy channel histograms");
577 Bool_t status = kTRUE;
578 AliITSdigitSPD *digitSPD = 0;
581 TBranch *digitsSPDBranch;
585 fITSLoader->LoadDigits("read");
587 // Create noisy channel histograms
588 AliITSPreprocessorSPD::CreateHistograms();
590 // Create an empty SPD digits array
591 TClonesArray *digitsArraySPD = AliITSPreprocessorSPD::CreateDigitsArray();
593 // Get number of events
594 UInt_t numberOfEvents = (fRunLoader->TreeE()) ? static_cast<UInt_t>(fRunLoader->TreeE()->GetEntries()) : 0;
596 // Make sure we don't try to analyze more data than there actually is
597 if (numberOfEvents < fMaximumNumberOfEvents)
599 fMaximumNumberOfEvents = numberOfEvents;
602 // Loop over all digits and put them in the corresponding histograms
603 if (numberOfEvents > 0)
605 for (UInt_t event = 0; event < fMaximumNumberOfEvents; event++)
607 // Get the current event
608 fRunLoader->GetEvent(event);
610 // Get the ITS digits tree
611 digitsTree = fITSLoader->TreeD();
613 // Disable all branches except the SPD branch to speed up the reading process
614 digitsTree->SetBranchStatus("ITSDigitsSPD",1);
615 digitsTree->SetBranchStatus("ITSDigitsSDD",0);
616 digitsTree->SetBranchStatus("ITSDigitsSSD",0);
618 // Reset the SPD digits array to make sure it is empty
619 digitsArraySPD->Clear();
621 // Get the SPD digits branch and set the address
622 digitsSPDBranch = digitsTree->GetBranch("ITSDigitsSPD");
624 digitsSPDBranch->SetAddress(&digitsArraySPD);
626 if (event%100 == 0) AliInfo(Form("Event #%d", event));
628 // Loop over all modules
629 UInt_t numberOfDigits = 0;
630 for (UInt_t module = 0; module < fNumberOfModules; module++)
632 // Get event data for current module
633 digitsTree->GetEvent(module);
635 // Get the number of entries
636 numberOfDigits = digitsArraySPD->GetEntries();
638 // Loop over all digits
639 for (UInt_t digit = 0; digit < numberOfDigits; digit++)
641 // Get the current digit
642 digitSPD = (AliITSdigitSPD*) digitsArraySPD->At(digit);
643 column = digitSPD->GetCoord1();
644 row = digitSPD->GetCoord2();
646 // Fill the digits histogram
647 ((TH2F*)(*fDigitsHistogram)[module])->Fill(column, row);
652 } // end numberOfEvents > 0
660 delete digitsArraySPD;
667 //__________________________________________________________________________
668 Bool_t AliITSPreprocessorSPD::FindDeadChannels(void)
670 // Locate dead channels
672 // Input : Filled hit histograms
673 // Output: TObjArray (fBadChannelsObjArray) with all identified dead channels
674 // Return: kTRUE if dead channels have been found
676 Bool_t status = kFALSE;
678 // Proceed only if properly initialized
683 // Initialize counters
684 fNumberOfBadChannels = new Int_t[fNumberOfModules];
685 for (UInt_t module = 0; module < fNumberOfModules; module++)
687 fNumberOfBadChannels[module] = 0;
690 // Loop over all modules (intentional modularization - DCS will occationally want to
691 // look for noisy channels in certain modules, but not all)
692 fIndex = 0; // Global bad channels array counter (must be reset here)
693 for (UInt_t module = 0; module < fNumberOfModules; module++)
695 status |= AliITSPreprocessorSPD::FindDeadChannelsInModule(module);
700 AliError("Dead channels can only be found in data taken with stand-alone VME system");
705 AliError("Not properly initialized");
712 //__________________________________________________________________________
713 Bool_t AliITSPreprocessorSPD::FindDeadChannelsInModule(UInt_t module)
715 // Locate dead channels
716 // This method assumes that the preprocessor is operator in VME mode.
717 // The algorithm is very simple. It assumes that the data was taken in
718 // a mode where all working SPD pixels should respond as being hit.
719 // fThreshold is used as the limit where everything below this value will
720 // be considered as "dead".
722 // Input : Filled hit histograms
723 // Output: TObjArray (fBadChannelsObjArray) with all identified noisy channels
724 // Return: kTRUE if dead channels have been found
726 // Store the index number for this module
727 fBadChannelsIndexArray[module] = fIndex++;
729 UInt_t xBin, numberOfXBins;
730 UInt_t yBin, numberOfYBins;
733 numberOfXBins = ((TH2F*)(*fDigitsHistogram)[module])->GetNbinsX();
734 numberOfYBins = ((TH2F*)(*fDigitsHistogram)[module])->GetNbinsY();
736 // Loop over all bins in this histogram
737 for (xBin = 1; xBin <= numberOfXBins; xBin++)
738 for (yBin = 1; yBin <= numberOfYBins; yBin++)
740 binContent = ((TH2F*)(*fDigitsHistogram)[module])->GetBinContent(xBin, yBin);
742 // Store this channel/bin if outside accepted limits
743 // A channel has to fire MORE times than the fThreshold value, or it will
744 // be considered as "dead"
745 if (binContent <= fThreshold)
747 // Store the dead channel in the array
748 // The channel object will be deleted in the destructor using the TObjArray Delete() method
749 // (The array will assume ownership of the channel object)
750 AliITSChannelSPD *channel = new AliITSChannelSPD((Int_t)(xBin - 1), (Int_t)(yBin - 1));
752 // Store the noisy channel in the array
753 fBadChannelsObjArray->Add(channel);
755 // Keep track of the number of bad channels in this module
756 fNumberOfBadChannels[module]++;
759 // Keep track of the highest module number
760 if (module > fHighestModuleNumber) fHighestModuleNumber = module;
762 //AliInfo(Form("New dead pixel in (m,c,r) = (%d,%d,%d)", module, xBin - 1, yBin - 1));
766 return (fNumberOfBadChannels[module] > 0);
770 //__________________________________________________________________________
771 Bool_t AliITSPreprocessorSPD::FindNoisyChannels(void)
773 // Locate noisy channels by searching through the digit histograms
774 // (There is one digit histogram per SPD module, i.e. a half-stave, 5 chips)
777 // Output: TObjArray (fBadChannelsObjArray) with all identified noisy channels
778 // Return: kTRUE if noisy channels have been found
780 Bool_t status = kFALSE;
782 // Proceed only if properly initialized
785 // (For testing purposes, noisy channels can be inserted here)
786 //SetFakeNoisyChannel(4,10,10);
787 // Initialize counters
788 fNumberOfBadChannels = new Int_t[fNumberOfModules];
789 for (UInt_t module = 0; module < fNumberOfModules; module++)
791 fNumberOfBadChannels[module] = 0;
794 // Scan through all the histogram bins and search for average filling deviations
795 if ((fSelectedAlgorithm == kOptimizedForRealData) || (fSelectedAlgorithm == kOptimizedForRealDataRMS))
797 // Noisy channel algorithm optimized for real data..........................
798 // Histograms can have any shape (both thresholds and quotients are used)
799 // This algorithm can be used to find noisy channels even if the data was
800 // taken with beam. All channels outside accepted limits (set by fThreshold
801 // and fThresholdRatio) will be marked as noisy
803 if (fSelectedAlgorithm == kOptimizedForRealData)
805 AliInfo("Searching for noisy channels (optimized for real data)");
809 AliInfo("Searching for noisy channels (optimized for real data, RMS version)");
812 // Loop over all modules (intentional modularization - DCS will occationally want to
813 // look for noisy channels in certain modules, but not all)
814 fIndex = 0; // Global bad channels array counter (must be reset here)
815 for (UInt_t module = 0; module < fNumberOfModules; module++)
817 status |= AliITSPreprocessorSPD::FindNoisyChannelsInModuleAlgo0(module);
822 // Noisy channel algorithm optimized for calibration data...........................
823 // Histograms will/should only contain noisy channels (only thresholds are used)
824 // Calibration histograms should have no background. The calibration data
825 // should have been taken without beam. All channels outside accepted limits
826 // (set by fThreshold) will be marked as noisy
828 AliInfo("Searching for noisy channels (optimized for calibration data)");
830 // Loop over all modules (intentional modularization - DCS will occationally want to
831 // look for noisy channels in certain modules, but not all)
832 fIndex = 0; // Global bad channels array counter (must be reset here)
833 for (UInt_t module = 0; module < fNumberOfModules; module++)
835 status |= AliITSPreprocessorSPD::FindNoisyChannelsInModuleAlgo1(module);
841 AliError("Not properly initialized");
848 //__________________________________________________________________________
849 Bool_t AliITSPreprocessorSPD::FindNoisyChannelsInModuleAlgo0(UInt_t module)
851 // Locate the noisy channels in a module (optimized for real data)
853 // For each channel in these histograms, the algorithm checks the average
854 // filling of the neighboring 3, 5 or 8 channels (depending on the location
855 // of the channel in question; corner, border or inside), or compares with the
856 // RMS of the neighbors. If the average is "much" lower, the channel will be
857 // considered to be noisy. The default noisy-to-normal fraction is stored in the
858 // fThresholdRatio varible. It can be set with the SetThresholdRatio method.
859 // The channel also has to have fired more than a certain minimum, fThreshold.
860 // It can be set with the SetThreshold method.
862 // To avoid difficulties with noisy channels that occur in pairs, the
863 // neighboring channel with largest number of fillings will be removed from
864 // the average calculation.
866 // NOTE: Since this method modifies the fBadChannelsObjArray and fBadChannelsIndexArray
867 // it is essential to initialize the fIndex counter before calling this module
868 // the first time. The bad channel data does not have to be ordered per module
869 // in the fBadChannelsObjArray, but the indices of where the data of a certain module
870 // starts has to be correct. A wrong fIndex can lead to segmentation violation
872 // Input : module number, filled digit histograms
873 // Output: TObjArray (fBadChannelsObjArray) with all identified noisy channels,
874 // Int_t[] (fBadChannelsIndexArray) with fBadChannelsObjArray module indices,
875 // number of noisy channels in this module (global variable fNumberOfBadChannels[module])
876 // Return: kTRUE if there are noisy channels in this module
878 // Store the index number for this module
879 fBadChannelsIndexArray[module] = fIndex++;
881 UInt_t xBin, numberOfXBins;
882 UInt_t yBin, numberOfYBins;
885 UInt_t numberOfNeighboringBins;
887 Float_t sumBinContent;
888 Float_t neighborBinContent;
889 Float_t maxBinContent;
890 Float_t averageBinContent;
893 numberOfXBins = (UInt_t)((TH2F*)(*fDigitsHistogram)[module])->GetNbinsX();
894 numberOfYBins = (UInt_t)((TH2F*)(*fDigitsHistogram)[module])->GetNbinsY();
896 // Loop over all bins in this histogram
897 for (xBin = 1; xBin <= numberOfXBins; xBin++)
898 for (yBin = 1; yBin <= numberOfYBins; yBin++)
900 numberOfNeighboringBins = 0;
901 averageBinContent = 0.;
903 binContent = ((TH2F*)(*fDigitsHistogram)[module])->GetBinContent(xBin, yBin);
906 // Calculate the average pixel level on the surrounding pixels
907 for (neighborXBin = xBin - 1; neighborXBin <= xBin + 1; neighborXBin++)
908 for (neighborYBin = yBin - 1; neighborYBin <= yBin + 1; neighborYBin++)
910 if ( !((neighborXBin == xBin) && (neighborYBin == yBin)) )
912 // Only update the number of neighboring bins when we are not on the border
913 if ((neighborXBin > 0) && (neighborXBin <= numberOfXBins+1) &&
914 (neighborYBin > 0) && (neighborYBin <= numberOfYBins+1))
917 ((TH2F*)(*fDigitsHistogram)[module])->GetBinContent(neighborXBin, neighborYBin);
919 if (fSelectedAlgorithm == kOptimizedForRealDataRMS)
922 sumBinContent += neighborBinContent*neighborBinContent;
927 sumBinContent += neighborBinContent;
930 if (neighborBinContent > maxBinContent) maxBinContent = neighborBinContent;
932 numberOfNeighboringBins++;
937 // Calculate the average. Remove the largest neighboring bin
938 // (Correction for potential clusters of noisy channels)
939 if (fSelectedAlgorithm == kOptimizedForRealDataRMS)
941 // Square the max bin content before removing it from the average calculation
942 maxBinContent *= maxBinContent;
944 averageBinContent = TMath::Sqrt((sumBinContent - maxBinContent)/(Float_t)(numberOfNeighboringBins - 1));
949 averageBinContent = (sumBinContent - maxBinContent)/(Float_t)(numberOfNeighboringBins - 1);
952 // Store this channel/bin if outside accepted limits
953 // The threshold ratio is the threshold for the current bin content divided by the
954 // average neighboring bin contents. The threshold bin content is the minimum number of
955 // times a channel has to have fired to be called noisy
956 ratio = (averageBinContent > 0) ? binContent/averageBinContent : 0.;
957 if ( ((ratio >= fThresholdRatio) || (ratio == 0.)) && (binContent >= fThreshold) )
959 // Store the noisy channel in the array
960 // The channel object will be deleted in the destructor using the TObjArray Delete() method
961 // (The array will assume ownership of the channel object)
962 AliITSChannelSPD *channel = new AliITSChannelSPD((Int_t)(xBin - 1), (Int_t)(yBin - 1));
964 // Store the noisy channel in the array
965 fBadChannelsObjArray->Add(channel);
967 // Keep track of the number of bad channels in this module
968 fNumberOfBadChannels[module]++;
971 // Keep track of the highest module number
972 if (module > fHighestModuleNumber) fHighestModuleNumber = module;
974 AliInfo(Form("New noisy pixel in (m,c,r) = (%d,%d,%d)", module, xBin - 1, yBin - 1));
975 AliInfo(Form("- Noisy pixel fired %d times, average neighborhood: %f",(Int_t)binContent,averageBinContent));
979 return (fNumberOfBadChannels[module] > 0);
983 //__________________________________________________________________________
984 Bool_t AliITSPreprocessorSPD::FindNoisyChannelsInModuleAlgo1(UInt_t module)
986 // Locate the noisy channels in a module (optimized for calibration data)
988 // This algorithm locates noisy channels by assuming original data was taken
989 // in calibration mode. This should be done without beam and will thus only
990 // contain data corresponding to background and noisy channels. The latter
991 // should be clearly visible in this data so this algorithm simply assumes
992 // that all histogram bins that are filled more than fThreshold times are
995 // NOTE: Since this method modifies the fBadChannelsObjArray and fBadChannelsIndexArray
996 // it is essential to initialize the fIndex counter before calling this module
997 // the first time. The bad channel data does not have to be ordered per module
998 // in the fBadChannelsObjArray, but the indices of where the data of a certain module
999 // starts has to be correct. A wrong fIndex can lead to segmentation violation
1001 // Input : module number, filled digit histograms
1002 // Output: TObjArray (fBadChannelsObjArray) with all identified noisy channels,
1003 // Int_t[] (fBadChannelsIndexArray) with fBadChannelsObjArray module indices,
1004 // number of noisy channels in this module (global variable fNumberOfBadChannels[module])
1005 // Return: kTRUE if there are noisy channels in this module
1007 // Store the index number for this module
1008 fBadChannelsIndexArray[module] = fIndex++;
1010 UInt_t xBin, numberOfXBins;
1011 UInt_t yBin, numberOfYBins;
1014 numberOfXBins = ((TH2F*)(*fDigitsHistogram)[module])->GetNbinsX();
1015 numberOfYBins = ((TH2F*)(*fDigitsHistogram)[module])->GetNbinsY();
1017 // Loop over all bins in this histogram
1018 for (xBin = 1; xBin <= numberOfXBins; xBin++)
1019 for (yBin = 1; yBin <= numberOfYBins; yBin++)
1021 binContent = ((TH2F*)(*fDigitsHistogram)[module])->GetBinContent(xBin, yBin);
1023 // Store this channel/bin if outside accepted limits
1024 // The threshold bin content is the minimum number of times a channel has to have
1025 // fired to be called noisy
1026 if (binContent >= fThreshold)
1028 // Store the noisy channel in the array
1029 // The channel object will be deleted in the destructor using the TObjArray Delete() method
1030 // (The array will assume ownership of the channel object)
1031 AliITSChannelSPD *channel = new AliITSChannelSPD((Int_t)(xBin - 1), (Int_t)(yBin - 1));
1033 // Store the noisy channel in the array
1034 fBadChannelsObjArray->Add(channel);
1036 // Keep track of the number of bad channels in this module
1037 fNumberOfBadChannels[module]++;
1040 // Keep track of the highest module number
1041 if (module > fHighestModuleNumber) fHighestModuleNumber = module;
1043 AliInfo(Form("New noisy pixel in (m,c,r) = (%d,%d,%d)", module, xBin - 1, yBin - 1));
1044 AliInfo(Form("- Noisy pixel fired %d times",(Int_t)binContent));
1048 return (fNumberOfBadChannels[module] > 0);
1052 //__________________________________________________________________________
1053 void AliITSPreprocessorSPD::PrintChannels(void)
1055 // Print all found bad channels to stdout
1057 // Input : fBadChannelsObjArray
1058 // Output: (dump to stdout)
1063 AliITSChannelSPD *channel = 0;
1065 // Print the bad channels stores in the array
1066 AliInfo("\nModule #\tColumn #\tRow #\n------------------------------------------------");
1067 for (UInt_t module = 0; module < fNumberOfModules; module++)
1070 while (j < fNumberOfBadChannels[module])
1072 channel = (AliITSChannelSPD *) fBadChannelsObjArray->At(i++);
1073 std::cout << module << "\t\t" << channel->GetColumn() << "\t\t" << channel->GetRow() << std::endl;
1075 // Go to next bad channel
1080 AliInfo(Form("%d bad channels were found", fBadChannelsObjArray->GetEntries()));
1084 //__________________________________________________________________________
1085 void AliITSPreprocessorSPD::MarkNoisyChannels(void)
1087 // WARNING: THIS METHOD DOESN'T WORK!!!
1089 // Mark all identified noisy channels
1091 // Input : List of noisy channels, original digits tree
1092 // Output: New digits tree containing SPD digits marked when noisy
1095 // The original digits tree (digitsTree) is cloned except for the SPD branch (ITSDigitSPD).
1096 // This branch is then redefined for each event and will contain all the original
1097 // information. All known noisy channels will be marked by using the TObject status bits
1098 // according to the following scheme. Dead channels are included for completeness. Note
1099 // that a dead channel will NEVER show up among digits..
1101 // Interpretation of digit status bits (LSB):
1102 // Dead channel Noisy channel | Integer
1103 // -----------------------------------------
1108 // meaning e.g. that a channel that is noisy will have the first bit set in its status bits
1110 // Do not continue unless we are processing DAQ data
1113 AliInfo("Marking bad channels");
1115 // Create the storage container that will be used to access the bad channels
1116 if (!fBadChannelsContainer)
1118 // Add the bad channels array to the storage container
1119 // (ownership is passed to the AliRunDataStorage object)
1120 fBadChannelsContainer = new AliITSBadChannelsSPD();
1122 // Convert the bad channels from TObjArray to Int_t[]
1123 AliITSPreprocessorSPD::ConvertObjToIntArray();
1125 // Store the arrays in the bad channels container object
1126 const Int_t kBadChannelsArraySize =
1127 2*fBadChannelsObjArray->GetEntries() + fNumberOfModules;
1128 fBadChannelsContainer->Put(fBadChannelsIntArray, kBadChannelsArraySize,
1129 fBadChannelsIndexArray, fNumberOfModules);
1132 // Create the bad channels helper object
1133 // (will be used to find a bad channel within a TObjArray)
1134 AliITSBadChannelsAuxSPD *aux = new AliITSBadChannelsAuxSPD();
1136 AliITSdigitSPD *digitSPD = 0;
1137 UInt_t numberOfDigits;
1139 Bool_t mark = kFALSE;
1141 TBranch *digitsBranch = 0;
1144 // Create an empty SPD digit array
1145 TObjArray *digitsArraySPD = new TObjArray();
1147 // Get the digits in update mode (we want to modify them if there are noisy channels)
1148 fITSLoader->UnloadDigits();
1149 fITSLoader->LoadDigits("update");
1151 // Get the number of events
1152 UInt_t numberOfEvents = (fRunLoader->TreeE()) ? static_cast<UInt_t>(fRunLoader->TreeE()->GetEntries()) : 0;
1154 // Loop over all events
1155 for (UInt_t event = 0; event < numberOfEvents; event++)
1157 if (event%100 == 0) AliInfo(Form("Event #%d", event));
1159 // Get the current event
1160 fRunLoader->GetEvent(event);
1162 // Get the ITS digits tree
1163 digitsTree = fITSLoader->TreeD();
1165 // Get SPD branch that will contain all digits with marked noisy channels
1166 digitsBranch = digitsTree->GetBranch("ITSDigitsSPD");
1167 digitsBranch->SetAddress(&digitsArraySPD);
1169 // Get the stored number of modules
1170 UInt_t numberOfModules = (Int_t)digitsTree->GetEntries();
1171 TObjArray **newDigitsArraySPD = new TObjArray*[numberOfModules];
1173 Int_t *digitNumber = new Int_t[numberOfModules];
1174 for (UInt_t m = 0; m < numberOfModules; m++)
1176 newDigitsArraySPD[m] = new TObjArray();
1180 AliInfo(Form("ent = %d", (Int_t)digitsTree->GetEntries()));
1182 // Reset the SPD digit arrays to make sure they are empty
1183 digitsArraySPD->Clear();
1185 // Get the SPD digits branch from the original digits tree and set the address
1186 digitsBranch = digitsTree->GetBranch("ITSDigitsSPD");
1187 digitsBranch->SetAddress(&digitsArraySPD);
1189 // Loop over all modules
1190 for (UInt_t module = 0; module < fNumberOfModules; module++)
1192 // Get event data for current module
1193 digitsTree->GetEvent(module);
1195 // Get the hits in the current module
1196 TObjArray *moduleObjArray = fBadChannelsContainer->CreateModuleObjArray(module);
1198 // Get the number of entries
1199 numberOfDigits = digitsArraySPD->GetEntries();
1201 // Loop over all digits and all channels
1202 for (UInt_t digit = 0; digit < numberOfDigits; digit++)
1204 // Get the current digit
1205 digitSPD = (AliITSdigitSPD*) digitsArraySPD->At(digit);
1206 newDigit[0] = digitSPD->GetCoord1(); // column
1207 newDigit[1] = digitSPD->GetCoord2(); // row
1208 newDigit[2] = digitSPD->GetSignal(); // signal
1210 // Check if this channel is noisy
1211 // (Compare with all stored channels in the bad channels array)
1212 if (aux->Find(digitSPD, moduleObjArray))
1214 // Set the mark flag and break the loop
1218 // Store this digit in the SPD digits array using a placement new operation
1219 new ((*newDigitsArraySPD[module])[digitNumber[module]]) AliITSdigitSPD(newDigit);
1221 // Mark it if noisy and store in the noisy channel array
1224 // Store this digit in the marked SPD digits array using a placement new operation
1225 //new ((*badChannels[m])[numberOfBadChannels[m]]) AliITSChannelSPD(newBadChannel);
1226 //new ((*newDigitsArraySPD[module])[digitNumber[module]]) AliITSdigitSPD(newDigit);
1228 // Mark the original channel as noisy
1229 ((*newDigitsArraySPD[module])[digitNumber[module]])->SetBit(kNoisyChannel);
1234 digitNumber[module]++;
1239 delete moduleObjArray;
1242 } // end module loop
1244 digitsBranch->Reset();
1245 digitsBranch->ResetAddress();
1248 delete digitsArraySPD;
1250 digitsTree->Reset();
1252 // WHY THIS RANGE?????????????????????????????????????????????????????????????????????
1253 for (UInt_t n = 0; n < event; n++)
1255 digitsTree->SetBranchAddress("ITSDigitsSPD", &newDigitsArraySPD[n]);
1259 digitsTree->AutoSave();
1262 for (UInt_t n = 0; n < event; n++)
1264 delete newDigitsArraySPD[n];
1266 delete [] newDigitsArraySPD;
1267 newDigitsArraySPD = 0;
1268 delete [] digitNumber;
1273 } // end loop over all events
1275 // Unload the digits
1276 fITSLoader->UnloadDigits();
1286 //__________________________________________________________________________
1287 Bool_t AliITSPreprocessorSPD::Store(AliCDBId &id, AliCDBMetaData *md)
1289 // Store the bad channels object in the calibration database
1290 // (See the corresponding run macro for further explanations)
1292 // Input : fBadChannelsObjArray (now containing all found bad channels), object meta data
1293 // Output: Database file containing the bad channels
1294 // Return: kTRUE if successful
1296 Bool_t status = kFALSE;
1298 AliInfo("Storing bad channels");
1300 // Add the bad channels array to the storage container
1301 // (ownership is passed to the AliRunDataStorage object)
1302 fBadChannelsContainer = new AliITSBadChannelsSPD();
1304 // Convert the bad channels from TObjArray to Int_t[]
1305 AliITSPreprocessorSPD::ConvertObjToIntArray();
1307 // Store the arrays in the bad channels container object
1308 const Int_t kBadChannelsArraySize =
1309 2*fBadChannelsObjArray->GetEntries() + fNumberOfModules;
1310 fBadChannelsContainer->Put(fBadChannelsIntArray, kBadChannelsArraySize,
1311 fBadChannelsIndexArray, fNumberOfModules);
1313 // Store the container
1314 if(!AliCDBManager::Instance()->IsDefaultStorageSet()) {
1315 //AliError("No storage set!");
1317 AliCDBManager::Instance()->SetDefaultStorage("local://Calib");
1320 if (AliCDBManager::Instance()->GetDefaultStorage()->Put(fBadChannelsContainer, id, md))
1322 AliInfo("Bad channels object stored in database");
1327 AliError("Failed to store object in database");
1334 //__________________________________________________________________________
1335 void AliITSPreprocessorSPD::ConvertObjToIntArray()
1337 // Convert the bad channel TObjArray to an Int_t array
1339 // Input : fBadChannelsObjArray (now containing all found bad channels)
1340 // Output: fBadChannelsIntArray
1344 // The TObjArray of this class (fBadChannelsObjArray) is converted to a sequential
1345 // Int_t array (fBadChannelsIntArray) in this method. For each module, the first
1346 // stored number is the number of bad channels in the current module. This is
1347 // followed by all the columns and rows of the bad channels:
1349 // badChannelsArray =
1350 // | N(m) | col0 | row0 | .. | colN(m) | N(m+1) | col0 | row0 | ...
1351 // . .......... module m ......... . .... module m+1 ......
1353 // The bad channels index array (fBadChannelsIndexArray) contains the indices of
1354 // the badChannelsArray, i.e. where the bad channels in certain module starts:
1356 // fBadChannelsObjArray =
1357 // | i0 | i1 | .. | iM | (where M = the number of SPD modules)
1359 // e.g. i1 corresponds to the index of the badChannelsArray where N(1) is stored,
1360 // i.e. the number of bad channels for module 1
1362 const Int_t kBadChannelsArraySize =
1363 2*fBadChannelsObjArray->GetEntries() + fNumberOfModules;
1364 fBadChannelsIntArray = new Int_t[kBadChannelsArraySize]; // Will be deleted in dtor
1365 AliITSChannelSPD *channel = 0;
1370 // Loop over all modules
1371 for (UInt_t module = 0; module < fNumberOfModules; module++)
1373 // Encode the number of bad channels of the current module
1374 fBadChannelsIntArray[k++] = fNumberOfBadChannels[module];
1376 // The columns and rows of the fBadChannelsObjArray will be stored sequentially
1377 // in the Int_t array
1379 while (j < fNumberOfBadChannels[module])
1381 channel = (AliITSChannelSPD *) fBadChannelsObjArray->At(i++);
1382 fBadChannelsIntArray[k++] = channel->GetColumn();
1383 fBadChannelsIntArray[k++] = channel->GetRow();
1385 // Go to next bad channel
1392 //__________________________________________________________________________
1393 Bool_t AliITSPreprocessorSPD::Store(AliCDBId& /*id*/, AliCDBMetaData* /*md*/, Int_t runNumber)
1395 // Store the bad channels object in the calibration database
1396 // (See the corresponding run macro for further explanations)
1398 // Input : fBadChannelsObjArray (now containing all found bad channels), object meta data
1399 // Output: Database file containing the bad channels
1400 // Return: kTRUE if successful
1402 Bool_t status = kFALSE;
1404 AliInfo("Storing bad channels");
1406 if(!AliCDBManager::Instance()->IsDefaultStorageSet()) {
1407 AliWarning("No storage set! Will use dummy one");
1408 AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
1412 AliCDBEntry *entrySPD = AliCDBManager::Instance()->Get("ITS/Calib/CalibSPD", runNumber);
1414 AliWarning("Calibration object retrieval failed! Dummy calibration will be used.");
1415 AliCDBStorage *origStorage = AliCDBManager::Instance()->GetDefaultStorage();
1416 AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
1418 entrySPD = AliCDBManager::Instance()->Get("ITS/Calib/CalibSPD", runNumber);
1419 AliCDBManager::Instance()->SetDefaultStorage(origStorage);
1422 TObjArray *respSPD = (TObjArray *)entrySPD->GetObject();
1425 AliWarning("Can not get calibration from calibration database !");
1431 AliITSChannelSPD *channel = 0;
1432 AliITSCalibrationSPD* res;
1433 for (Int_t module=0; module<respSPD->GetEntries(); module++) {
1434 res = (AliITSCalibrationSPD*) respSPD->At(module);
1436 while (j < fNumberOfBadChannels[module]) {
1437 channel = (AliITSChannelSPD *) fBadChannelsObjArray->At(i++);
1438 res->AddDead(channel->GetColumn(),channel->GetRow());
1439 // Go to next bad channel
1445 AliCDBManager::Instance()->Put(entrySPD);
1446 entrySPD->SetObject(NULL);
1447 entrySPD->SetOwner(kTRUE);