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.3 2006/04/12 08:32:31 hristov
19 New SPD simulation (Massimo):
20 - speeding up of the diffusion code (Bjorn)
21 - creation of CDB file with the dead channels, implementation
22 of the CDB reading, check of the code (Henrik, Giuseppe, Domenico)
23 - final tuning of the diffusion model parameters (Romualdo)
25 Revision 1.2 2006/02/03 11:31:18 masera
26 Calibration framework improved (E. Crescio)
28 Revision 1.1 2005/10/11 12:31:50 masera
29 Preprocessor classes for SPD (Paul Nilsson)
33 ///////////////////////////////////////////////////////////////////////////
34 // AliITSPreprocessorSPD implementation by P. Nilsson 2005
35 // MAIN AUTHOR/CONTACT: Paul.Nilsson@cern.ch
37 // The purpose of this class is to provide algorithms for identification
38 // of "bad channels" such as dead channels and noisy channels in the SPD
40 // Examples on how to use this class can be found in the macros:
42 // .findNoisyChannels.C (Locate and store noisy channels in the CDB)
43 // .readNoisyChannels.C (Read noisy channels from the CDB)
44 // .findDeadChannels.C (Locate and store dead channels in the CDB)
45 // .readDeadChannels.C (Read dead channels from the CDB)
47 // Modified by D. Elia, H. Tydesjo
48 // March 2006: Mixed up coordinates, bug fixed
50 ///////////////////////////////////////////////////////////////////////////
53 #include "AliITSPreprocessorSPD.h"
54 #include "AliCDBEntry.h"
55 #include "AliITSCalibrationSPD.h"
56 ClassImp(AliITSPreprocessorSPD)
59 //__________________________________________________________________________
60 AliITSPreprocessorSPD::AliITSPreprocessorSPD(void):
65 fMaximumNumberOfEvents(1000000),
66 fHighestModuleNumber(0),
67 fSelectedAlgorithm(kOptimizedForRealData),
68 fGeometryMode(kALICEGeometry),
69 fNumberOfBadChannels(0),
73 fBadChannelsObjArray(0),
74 fBadChannelsIntArray(0),
75 fBadChannelsIndexArray(0),
76 fBadChannelsContainer(0)
78 // Default constructor for the SPD preprocessor
80 // Initialization has to be done by hand using Set* methods and the Open method
84 // Return: <empty/uninitialized AliITSPreprocessorSPD object>
88 //__________________________________________________________________________
89 AliITSPreprocessorSPD::AliITSPreprocessorSPD(const char *fileName, const char *mode,
90 const char *fileNameg, const Int_t maxNumberOfEvents):
95 // Standard constructor for the SPD preprocessor
97 // Input : Name of digit file
99 // Return: Initialized SPD preprocessor object
102 AliITSPreprocessorSPD::Init();
103 AliITSPreprocessorSPD::SetMaximumNumberOfEvents(maxNumberOfEvents);
105 // Open and read the galice and digit files
106 if (!AliITSPreprocessorSPD::Open(fileName, mode, fileNameg))
108 AliError("Failed to open file");
114 //__________________________________________________________________________
115 AliITSPreprocessorSPD::AliITSPreprocessorSPD(const AliITSPreprocessorSPD &prep) :
118 // Default copy constructor
119 // Notice that only pointer addresses are copied!
120 // Memory allocations of new objects are not done.
122 fITSLoader = prep.fITSLoader;
123 fRunLoader = prep.fRunLoader;
124 fThresholdRatio = prep.fThresholdRatio;
125 fThreshold = prep.fThreshold;
126 fMaximumNumberOfEvents = prep.fMaximumNumberOfEvents;
127 fHighestModuleNumber = prep.fHighestModuleNumber;
128 fSelectedAlgorithm = prep.fSelectedAlgorithm;
129 fGeometryMode = prep.fGeometryMode;
130 fNumberOfBadChannels = prep.fNumberOfBadChannels;
132 fVMEMode = prep.fVMEMode;
133 fDigitsHistogram = prep.fDigitsHistogram;
134 fBadChannelsObjArray = prep.fBadChannelsObjArray;
135 fBadChannelsIntArray = prep.fBadChannelsIntArray;
136 fBadChannelsIndexArray = prep.fBadChannelsIndexArray;
137 fBadChannelsContainer = prep.fBadChannelsContainer;
141 //__________________________________________________________________________
142 AliITSPreprocessorSPD& AliITSPreprocessorSPD::operator=(const AliITSPreprocessorSPD &prep)
144 // Assignment operator
146 AliError("Not implemented");
148 if (this != &prep) { } // Do not delete this line
154 //__________________________________________________________________________
155 AliITSPreprocessorSPD::~AliITSPreprocessorSPD(void)
157 // Destructor for the SPD preprocessor
159 if (fDigitsHistogram)
161 // try fDigitsHistogram->Delete(); if the following crashes
162 for (UInt_t module = 0; module < fNumberOfModules; module++)
164 (*fDigitsHistogram)[module]->Delete();
166 delete fDigitsHistogram;
167 fDigitsHistogram = 0;
170 if (fNumberOfBadChannels)
172 delete [] fNumberOfBadChannels;
173 fNumberOfBadChannels = 0;
176 if (fBadChannelsIntArray)
178 delete [] fBadChannelsIntArray;
179 fBadChannelsIntArray = 0;
182 if (fBadChannelsIndexArray)
184 delete [] fBadChannelsIndexArray;
185 fBadChannelsIndexArray = 0;
188 if (fBadChannelsObjArray)
190 fBadChannelsObjArray->Delete();
191 delete fBadChannelsObjArray;
192 fBadChannelsObjArray = 0;
198 delete fBadChannelsContainer;
199 fBadChannelsContainer = 0;
203 //__________________________________________________________________________
204 void AliITSPreprocessorSPD::Init(void)
206 // Initialization of the SPD preprocessor
209 // Output: Several logistical variables
212 // Default maximum number of events per histogram
213 fMaximumNumberOfEvents = 1000000;
215 // Default noisy channel removal algorithm
216 fSelectedAlgorithm = kOptimizedForRealData;
218 // Default noisy channel threshold and threshold ratio
219 // (threshold for current bin content divided by the average neighboring bin contents)
221 fThresholdRatio = 5.;
223 // Set default geometry mode (ALICE geometry). This also sets fNumberOfModules
224 AliITSPreprocessorSPD::SetGeometryMode(kALICEGeometry);
226 // The highest module number with found bad channels
227 fHighestModuleNumber = 0;
229 // Assume input is not from a VME file
232 // Initialization is complete
237 //__________________________________________________________________________
238 void AliITSPreprocessorSPD::SetGeometryMode(UInt_t mode)
240 // Set the geometry mode
242 // Input : Geometry mode (either kTestBeamGeometry or kALICEGeometry)
246 fGeometryMode = mode;
248 // In case of an input VME file, the number of modules has already been fixed.
249 // Do not try to change it
252 if (mode == kALICEGeometry)
254 fNumberOfModules = kNumberOfSPDModules;
256 else if (mode == kTestBeamGeometry)
258 fNumberOfModules = kNumberOfTestBeamSPDModules;
262 AliError("Unknown geometry mode, defaults to ALICE geometry");
263 fNumberOfModules = kNumberOfSPDModules;
269 //__________________________________________________________________________
270 void AliITSPreprocessorSPD::SetFakeNoisyChannel(Int_t module, Int_t column, Int_t row)
272 // Introduce a fake noisy channel in the hit histograms
274 // Input : Module, column and row numbers
275 // Output: Updated hit histograms
278 if ((UInt_t)module < fNumberOfModules)
280 ((TH2F*)(*fDigitsHistogram)[module])->Fill(column, row, 1000000);
284 AliError("No such module number");
289 //__________________________________________________________________________
290 Bool_t AliITSPreprocessorSPD::Open(const char *fileName, const char *mode, const char *fileNameg)
294 // Input : Name and mode of ITS digit file, name of galice file
296 // Return: kTRUE if loaders succeed
298 Bool_t status = kFALSE;
299 Bool_t status0 = kFALSE;
300 Bool_t status1 = kFALSE;
301 Bool_t status2 = kFALSE;
302 Bool_t status3 = kFALSE;
304 // Only proceed if initialization has been done
307 TString m = (TString)mode;
308 if (m == "daq" || m == "DAQ")
310 // Open the data file and get the run loader
311 fRunLoader = AliRunLoader::Open(fileNameg);
314 // Get the gAlice object
315 status0 = AliITSPreprocessorSPD::GetgAlice();
317 // Get the ITS digits
318 if (status0) status1 = AliITSPreprocessorSPD::GetITSDigits(fileName);
320 // Create the test beam object
321 if (status1) status2 = AliITSPreprocessorSPD::CreateGeometryObj();
323 // Fill histograms with DAQ digits
324 if (status2) status3 = AliITSPreprocessorSPD::FillHistograms();
326 status = status0 & status1 & status2 & status3;
330 AliError("Failed to get the run loader");
333 else if (m == "vme" || m == "VME")
335 // Open the VME file that contains histograms with fired channels as read by the stand-alone VME system
336 TFile *vmeFile = TFile::Open(fileName);
338 if (!vmeFile->IsOpen())
340 AliError("Could not open VME input file");
344 // Get the histograms from the VME file that contains all fired channels
345 status0 = AliITSPreprocessorSPD::GetVMEHistograms(vmeFile);
347 // Create the test beam object
348 if (status0) status1 = AliITSPreprocessorSPD::CreateGeometryObj();
351 // This boolean will be used to override any attempts of changing the number of modules by the user
352 // with the SetGeometryMode method. For VME files the number of modules will entirely be determined by
353 // the input VME root file, i.e. by the number of histograms in this file
356 status = status0 & status1;
360 AliError("Unknown filetype - assuming DAQ file");
363 // At this stage, the final number of modules will be known (safe to define arrays)
364 // In case data is read from a VME root file, it will not be known before
367 // Create the arrays for bookkeeping and storing the noisy channels
368 if (!fBadChannelsObjArray)
370 fBadChannelsObjArray = new TObjArray();
372 if (!fBadChannelsIndexArray)
374 // This array will contain the begin and end indices for each module, i.e. where to begin
375 // and stop reading the fBadChannelsObjArray for a certain module.
376 // The last position of the array will contain the end index of the last module
377 fBadChannelsIndexArray = new Int_t[fNumberOfModules + 1];
378 for (UInt_t module = 0; module < fNumberOfModules + 1; module++)
380 fBadChannelsIndexArray[module] = 0;
387 AliError("SPD preprocessor not initialized. Can't load digits");
394 //__________________________________________________________________________
395 Bool_t AliITSPreprocessorSPD::GetgAlice(void)
397 // Get the gAlice object
400 // Output: gAlice object
401 // Return: kTRUE if the gAlice object was found
403 Bool_t status = kTRUE;
406 fRunLoader->LoadgAlice();
407 if (!fRunLoader->GetAliRun())
409 AliError("gAlice not found on file. Aborting.");
417 //__________________________________________________________________________
418 Bool_t AliITSPreprocessorSPD::GetVMEHistograms(TFile *vmeFile)
420 // Get pre-filled digit histograms from input VME file
422 // Input : pointer to VME file
424 // Return: kTRUE if histograms are found on file
426 Bool_t status = kFALSE;
428 // Get the file directory
429 TDirectory *dir = (TDirectory *)vmeFile;
431 // Get the number of keys (histograms in this case corresponding to modules/ladders)
432 fNumberOfModules = dir->GetNkeys();
433 if ((fNumberOfModules > 0) && (fNumberOfModules <= kNumberOfSPDModules))
437 // Create bad channel histograms
438 fDigitsHistogram = new TObjArray(fNumberOfModules);
440 // Create a key iterator
441 TIter nextkey(dir->GetListOfKeys());
444 // Loop over all objects, read them in to memory one by one
446 while ( (key = (TKey *)nextkey()) )
448 (*fDigitsHistogram)[module++] = (TH2F *)key->ReadObj();
451 if (module > fNumberOfModules)
453 AliError("Wrong number of keys in VME file");
461 AliError("Wrong number of histograms in VME file");
466 //__________________________________________________________________________
467 Bool_t AliITSPreprocessorSPD::GetITSDigits(const char *fileName)
469 // Get the ITS digits
471 // Input : name of digits file
472 // Output: fITSLoader object, ITS digits
473 // Return: kTRUE if loader succeed
475 Bool_t status = kTRUE;
477 // Load the gAlice and the header
478 fRunLoader->LoadgAlice();
479 fRunLoader->LoadHeader();
481 // Get the ITS loader
482 fITSLoader = (AliITSLoader*) fRunLoader->GetLoader("ITSLoader");
485 AliError("ITS loader not found");
490 // Open the digits file
491 fITSLoader->SetDigitsFileName(fileName);
498 //__________________________________________________________________________
499 TClonesArray* AliITSPreprocessorSPD::CreateDigitsArray(void) const
501 // Creation of the digits array
505 // Return: Pointer to the SPD digits array
507 // Create an array for 5 chips with 8192 channels each
508 TClonesArray *digitsArray = new TClonesArray("AliITSdigitSPD", kNumberOfChannels);
514 //__________________________________________________________________________
515 Bool_t AliITSPreprocessorSPD::CreateGeometryObj(void)
517 // Creation of the geometry object
519 // This object is used to get the number of SPD half-staves
522 // Output: fGeometryObj
523 // Return: kTRUE if fGeometryObj has been created
525 Bool_t status = true;
527 // Create geometry object
528 // fGeometryObj = new ...
529 // if (!fGeometryObj) status = kFALSE;
531 // Get SPD parameters
532 // fNumberOfColumns = fGeometryObject->GetNumberOfColumns();
533 // fNumberOfRows = fGeometryObject->GetNumberOfRows();
535 fNumberOfColumns = kNumberOfColumns;
536 fNumberOfRows = kNumberOfRows;
542 //__________________________________________________________________________
543 void AliITSPreprocessorSPD::CreateHistograms(void)
545 // Create the noisy channel histograms
548 // Output: Noisy channel histograms
552 char n[4]; // For storing the module number (maximum module number is 240)
554 // Create noisy channel histograms
555 fDigitsHistogram = new TObjArray(fNumberOfModules);
557 for (UInt_t i = 0; i < fNumberOfModules; i++)
560 moduleName = "module_";
562 moduleName.Append(n);
564 (*fDigitsHistogram)[i] = new TH2F(moduleName,"Digits",
565 fNumberOfColumns,-.5,(1.*fNumberOfColumns - .5),
566 fNumberOfRows,-.5,(1.*fNumberOfRows - .5));
571 //__________________________________________________________________________
572 Bool_t AliITSPreprocessorSPD::FillHistograms(void)
574 // Fill the histograms with digits (hit positions of unclustered hits)
576 // (There is one digit histogram per SPD module, i.e. a half-stave, 10 chips)
578 // Input : No arguments (Empty digit histograms)
579 // Output: Filled digit histograms
580 // Return: kTRUE if histograms are filled with digits, kFALSE otherwise
582 AliInfo("Filling noisy channel histograms");
584 Bool_t status = kTRUE;
585 AliITSdigitSPD *digitSPD = 0;
588 TBranch *digitsSPDBranch;
592 fITSLoader->LoadDigits("read");
594 // Create noisy channel histograms
595 AliITSPreprocessorSPD::CreateHistograms();
597 // Create an empty SPD digits array
598 TClonesArray *digitsArraySPD = AliITSPreprocessorSPD::CreateDigitsArray();
600 // Get number of events
601 UInt_t numberOfEvents = (fRunLoader->TreeE()) ? static_cast<UInt_t>(fRunLoader->TreeE()->GetEntries()) : 0;
603 // Make sure we don't try to analyze more data than there actually is
604 if (numberOfEvents < fMaximumNumberOfEvents)
606 fMaximumNumberOfEvents = numberOfEvents;
609 // Loop over all digits and put them in the corresponding histograms
610 if (numberOfEvents > 0)
612 for (UInt_t event = 0; event < fMaximumNumberOfEvents; event++)
614 // Get the current event
615 fRunLoader->GetEvent(event);
617 // Get the ITS digits tree
618 digitsTree = fITSLoader->TreeD();
620 // Disable all branches except the SPD branch to speed up the reading process
621 digitsTree->SetBranchStatus("ITSDigitsSPD",1);
622 digitsTree->SetBranchStatus("ITSDigitsSDD",0);
623 digitsTree->SetBranchStatus("ITSDigitsSSD",0);
625 // Reset the SPD digits array to make sure it is empty
626 digitsArraySPD->Clear();
628 // Get the SPD digits branch and set the address
629 digitsSPDBranch = digitsTree->GetBranch("ITSDigitsSPD");
631 digitsSPDBranch->SetAddress(&digitsArraySPD);
633 if (event%100 == 0) AliInfo(Form("Event #%d", event));
635 // Loop over all modules
636 UInt_t numberOfDigits = 0;
637 for (UInt_t module = 0; module < fNumberOfModules; module++)
639 // Get event data for current module
640 digitsTree->GetEvent(module);
642 // Get the number of entries
643 numberOfDigits = digitsArraySPD->GetEntries();
645 // Loop over all digits
646 for (UInt_t digit = 0; digit < numberOfDigits; digit++)
648 // Get the current digit
649 digitSPD = (AliITSdigitSPD*) digitsArraySPD->At(digit);
650 column = digitSPD->GetCoord1();
651 row = digitSPD->GetCoord2();
653 // Fill the digits histogram
654 ((TH2F*)(*fDigitsHistogram)[module])->Fill(column, row);
659 } // end numberOfEvents > 0
667 delete digitsArraySPD;
674 //__________________________________________________________________________
675 Bool_t AliITSPreprocessorSPD::FindDeadChannels(void)
677 // Locate dead channels
679 // Input : Filled hit histograms
680 // Output: TObjArray (fBadChannelsObjArray) with all identified dead channels
681 // Return: kTRUE if dead channels have been found
683 Bool_t status = kFALSE;
685 // Proceed only if properly initialized
690 // Initialize counters
691 fNumberOfBadChannels = new Int_t[fNumberOfModules];
692 for (UInt_t module = 0; module < fNumberOfModules; module++)
694 fNumberOfBadChannels[module] = 0;
697 // Loop over all modules (intentional modularization - DCS will occationally want to
698 // look for noisy channels in certain modules, but not all)
699 fIndex = 0; // Global bad channels array counter (must be reset here)
700 for (UInt_t module = 0; module < fNumberOfModules; module++)
702 status |= AliITSPreprocessorSPD::FindDeadChannelsInModule(module);
707 AliError("Dead channels can only be found in data taken with stand-alone VME system");
712 AliError("Not properly initialized");
719 //__________________________________________________________________________
720 Bool_t AliITSPreprocessorSPD::FindDeadChannelsInModule(UInt_t module)
722 // Locate dead channels
723 // This method assumes that the preprocessor is operator in VME mode.
724 // The algorithm is very simple. It assumes that the data was taken in
725 // a mode where all working SPD pixels should respond as being hit.
726 // fThreshold is used as the limit where everything below this value will
727 // be considered as "dead".
729 // Input : Filled hit histograms
730 // Output: TObjArray (fBadChannelsObjArray) with all identified noisy channels
731 // Return: kTRUE if dead channels have been found
733 // Store the index number for this module
734 fBadChannelsIndexArray[module] = fIndex++;
736 UInt_t xBin, numberOfXBins;
737 UInt_t yBin, numberOfYBins;
740 numberOfXBins = ((TH2F*)(*fDigitsHistogram)[module])->GetNbinsX();
741 numberOfYBins = ((TH2F*)(*fDigitsHistogram)[module])->GetNbinsY();
743 // Loop over all bins in this histogram
744 for (xBin = 1; xBin <= numberOfXBins; xBin++)
745 for (yBin = 1; yBin <= numberOfYBins; yBin++)
747 binContent = ((TH2F*)(*fDigitsHistogram)[module])->GetBinContent(xBin, yBin);
749 // Store this channel/bin if outside accepted limits
750 // A channel has to fire MORE times than the fThreshold value, or it will
751 // be considered as "dead"
752 if (binContent <= fThreshold)
754 // Store the dead channel in the array
755 // The channel object will be deleted in the destructor using the TObjArray Delete() method
756 // (The array will assume ownership of the channel object)
757 AliITSChannelSPD *channel = new AliITSChannelSPD((Int_t)(xBin - 1), (Int_t)(yBin - 1));
759 // Store the noisy channel in the array
760 fBadChannelsObjArray->Add(channel);
762 // Keep track of the number of bad channels in this module
763 fNumberOfBadChannels[module]++;
766 // Keep track of the highest module number
767 if (module > fHighestModuleNumber) fHighestModuleNumber = module;
769 //AliInfo(Form("New dead pixel in (m,c,r) = (%d,%d,%d)", module, xBin - 1, yBin - 1));
773 return (fNumberOfBadChannels[module] > 0);
777 //__________________________________________________________________________
778 Bool_t AliITSPreprocessorSPD::FindNoisyChannels(void)
780 // Locate noisy channels by searching through the digit histograms
781 // (There is one digit histogram per SPD module, i.e. a half-stave, 5 chips)
784 // Output: TObjArray (fBadChannelsObjArray) with all identified noisy channels
785 // Return: kTRUE if noisy channels have been found
787 Bool_t status = kFALSE;
789 // Proceed only if properly initialized
792 // (For testing purposes, noisy channels can be inserted here)
793 //SetFakeNoisyChannel(4,10,10);
794 // Initialize counters
795 fNumberOfBadChannels = new Int_t[fNumberOfModules];
796 for (UInt_t module = 0; module < fNumberOfModules; module++)
798 fNumberOfBadChannels[module] = 0;
801 // Scan through all the histogram bins and search for average filling deviations
802 if ((fSelectedAlgorithm == kOptimizedForRealData) || (fSelectedAlgorithm == kOptimizedForRealDataRMS))
804 // Noisy channel algorithm optimized for real data..........................
805 // Histograms can have any shape (both thresholds and quotients are used)
806 // This algorithm can be used to find noisy channels even if the data was
807 // taken with beam. All channels outside accepted limits (set by fThreshold
808 // and fThresholdRatio) will be marked as noisy
810 if (fSelectedAlgorithm == kOptimizedForRealData)
812 AliInfo("Searching for noisy channels (optimized for real data)");
816 AliInfo("Searching for noisy channels (optimized for real data, RMS version)");
819 // Loop over all modules (intentional modularization - DCS will occationally want to
820 // look for noisy channels in certain modules, but not all)
821 fIndex = 0; // Global bad channels array counter (must be reset here)
822 for (UInt_t module = 0; module < fNumberOfModules; module++)
824 status |= AliITSPreprocessorSPD::FindNoisyChannelsInModuleAlgo0(module);
829 // Noisy channel algorithm optimized for calibration data...........................
830 // Histograms will/should only contain noisy channels (only thresholds are used)
831 // Calibration histograms should have no background. The calibration data
832 // should have been taken without beam. All channels outside accepted limits
833 // (set by fThreshold) will be marked as noisy
835 AliInfo("Searching for noisy channels (optimized for calibration data)");
837 // Loop over all modules (intentional modularization - DCS will occationally want to
838 // look for noisy channels in certain modules, but not all)
839 fIndex = 0; // Global bad channels array counter (must be reset here)
840 for (UInt_t module = 0; module < fNumberOfModules; module++)
842 status |= AliITSPreprocessorSPD::FindNoisyChannelsInModuleAlgo1(module);
848 AliError("Not properly initialized");
855 //__________________________________________________________________________
856 Bool_t AliITSPreprocessorSPD::FindNoisyChannelsInModuleAlgo0(UInt_t module)
858 // Locate the noisy channels in a module (optimized for real data)
860 // For each channel in these histograms, the algorithm checks the average
861 // filling of the neighboring 3, 5 or 8 channels (depending on the location
862 // of the channel in question; corner, border or inside), or compares with the
863 // RMS of the neighbors. If the average is "much" lower, the channel will be
864 // considered to be noisy. The default noisy-to-normal fraction is stored in the
865 // fThresholdRatio varible. It can be set with the SetThresholdRatio method.
866 // The channel also has to have fired more than a certain minimum, fThreshold.
867 // It can be set with the SetThreshold method.
869 // To avoid difficulties with noisy channels that occur in pairs, the
870 // neighboring channel with largest number of fillings will be removed from
871 // the average calculation.
873 // NOTE: Since this method modifies the fBadChannelsObjArray and fBadChannelsIndexArray
874 // it is essential to initialize the fIndex counter before calling this module
875 // the first time. The bad channel data does not have to be ordered per module
876 // in the fBadChannelsObjArray, but the indices of where the data of a certain module
877 // starts has to be correct. A wrong fIndex can lead to segmentation violation
879 // Input : module number, filled digit histograms
880 // Output: TObjArray (fBadChannelsObjArray) with all identified noisy channels,
881 // Int_t[] (fBadChannelsIndexArray) with fBadChannelsObjArray module indices,
882 // number of noisy channels in this module (global variable fNumberOfBadChannels[module])
883 // Return: kTRUE if there are noisy channels in this module
885 // Store the index number for this module
886 fBadChannelsIndexArray[module] = fIndex++;
888 UInt_t xBin, numberOfXBins;
889 UInt_t yBin, numberOfYBins;
892 UInt_t numberOfNeighboringBins;
894 Float_t sumBinContent;
895 Float_t neighborBinContent;
896 Float_t maxBinContent;
897 Float_t averageBinContent;
900 numberOfXBins = (UInt_t)((TH2F*)(*fDigitsHistogram)[module])->GetNbinsX();
901 numberOfYBins = (UInt_t)((TH2F*)(*fDigitsHistogram)[module])->GetNbinsY();
903 // Loop over all bins in this histogram
904 for (xBin = 1; xBin <= numberOfXBins; xBin++)
905 for (yBin = 1; yBin <= numberOfYBins; yBin++)
907 numberOfNeighboringBins = 0;
908 averageBinContent = 0.;
910 binContent = ((TH2F*)(*fDigitsHistogram)[module])->GetBinContent(xBin, yBin);
913 // Calculate the average pixel level on the surrounding pixels
914 for (neighborXBin = xBin - 1; neighborXBin <= xBin + 1; neighborXBin++)
915 for (neighborYBin = yBin - 1; neighborYBin <= yBin + 1; neighborYBin++)
917 if ( !((neighborXBin == xBin) && (neighborYBin == yBin)) )
919 // Only update the number of neighboring bins when we are not on the border
920 if ((neighborXBin > 0) && (neighborXBin <= numberOfXBins+1) &&
921 (neighborYBin > 0) && (neighborYBin <= numberOfYBins+1))
924 ((TH2F*)(*fDigitsHistogram)[module])->GetBinContent(neighborXBin, neighborYBin);
926 if (fSelectedAlgorithm == kOptimizedForRealDataRMS)
929 sumBinContent += neighborBinContent*neighborBinContent;
934 sumBinContent += neighborBinContent;
937 if (neighborBinContent > maxBinContent) maxBinContent = neighborBinContent;
939 numberOfNeighboringBins++;
944 // Calculate the average. Remove the largest neighboring bin
945 // (Correction for potential clusters of noisy channels)
946 if (fSelectedAlgorithm == kOptimizedForRealDataRMS)
948 // Square the max bin content before removing it from the average calculation
949 maxBinContent *= maxBinContent;
951 averageBinContent = TMath::Sqrt((sumBinContent - maxBinContent)/(Float_t)(numberOfNeighboringBins - 1));
956 averageBinContent = (sumBinContent - maxBinContent)/(Float_t)(numberOfNeighboringBins - 1);
959 // Store this channel/bin if outside accepted limits
960 // The threshold ratio is the threshold for the current bin content divided by the
961 // average neighboring bin contents. The threshold bin content is the minimum number of
962 // times a channel has to have fired to be called noisy
963 ratio = (averageBinContent > 0) ? binContent/averageBinContent : 0.;
964 if ( ((ratio >= fThresholdRatio) || (ratio == 0.)) && (binContent >= fThreshold) )
966 // Store the noisy channel in the array
967 // The channel object will be deleted in the destructor using the TObjArray Delete() method
968 // (The array will assume ownership of the channel object)
969 AliITSChannelSPD *channel = new AliITSChannelSPD((Int_t)(xBin - 1), (Int_t)(yBin - 1));
971 // Store the noisy channel in the array
972 fBadChannelsObjArray->Add(channel);
974 // Keep track of the number of bad channels in this module
975 fNumberOfBadChannels[module]++;
978 // Keep track of the highest module number
979 if (module > fHighestModuleNumber) fHighestModuleNumber = module;
981 AliInfo(Form("New noisy pixel in (m,c,r) = (%d,%d,%d)", module, xBin - 1, yBin - 1));
982 AliInfo(Form("- Noisy pixel fired %d times, average neighborhood: %f",(Int_t)binContent,averageBinContent));
986 return (fNumberOfBadChannels[module] > 0);
990 //__________________________________________________________________________
991 Bool_t AliITSPreprocessorSPD::FindNoisyChannelsInModuleAlgo1(UInt_t module)
993 // Locate the noisy channels in a module (optimized for calibration data)
995 // This algorithm locates noisy channels by assuming original data was taken
996 // in calibration mode. This should be done without beam and will thus only
997 // contain data corresponding to background and noisy channels. The latter
998 // should be clearly visible in this data so this algorithm simply assumes
999 // that all histogram bins that are filled more than fThreshold times are
1002 // NOTE: Since this method modifies the fBadChannelsObjArray and fBadChannelsIndexArray
1003 // it is essential to initialize the fIndex counter before calling this module
1004 // the first time. The bad channel data does not have to be ordered per module
1005 // in the fBadChannelsObjArray, but the indices of where the data of a certain module
1006 // starts has to be correct. A wrong fIndex can lead to segmentation violation
1008 // Input : module number, filled digit histograms
1009 // Output: TObjArray (fBadChannelsObjArray) with all identified noisy channels,
1010 // Int_t[] (fBadChannelsIndexArray) with fBadChannelsObjArray module indices,
1011 // number of noisy channels in this module (global variable fNumberOfBadChannels[module])
1012 // Return: kTRUE if there are noisy channels in this module
1014 // Store the index number for this module
1015 fBadChannelsIndexArray[module] = fIndex++;
1017 UInt_t xBin, numberOfXBins;
1018 UInt_t yBin, numberOfYBins;
1021 numberOfXBins = ((TH2F*)(*fDigitsHistogram)[module])->GetNbinsX();
1022 numberOfYBins = ((TH2F*)(*fDigitsHistogram)[module])->GetNbinsY();
1024 // Loop over all bins in this histogram
1025 for (xBin = 1; xBin <= numberOfXBins; xBin++)
1026 for (yBin = 1; yBin <= numberOfYBins; yBin++)
1028 binContent = ((TH2F*)(*fDigitsHistogram)[module])->GetBinContent(xBin, yBin);
1030 // Store this channel/bin if outside accepted limits
1031 // The threshold bin content is the minimum number of times a channel has to have
1032 // fired to be called noisy
1033 if (binContent >= fThreshold)
1035 // Store the noisy channel in the array
1036 // The channel object will be deleted in the destructor using the TObjArray Delete() method
1037 // (The array will assume ownership of the channel object)
1038 AliITSChannelSPD *channel = new AliITSChannelSPD((Int_t)(xBin - 1), (Int_t)(yBin - 1));
1040 // Store the noisy channel in the array
1041 fBadChannelsObjArray->Add(channel);
1043 // Keep track of the number of bad channels in this module
1044 fNumberOfBadChannels[module]++;
1047 // Keep track of the highest module number
1048 if (module > fHighestModuleNumber) fHighestModuleNumber = module;
1050 AliInfo(Form("New noisy pixel in (m,c,r) = (%d,%d,%d)", module, xBin - 1, yBin - 1));
1051 AliInfo(Form("- Noisy pixel fired %d times",(Int_t)binContent));
1055 return (fNumberOfBadChannels[module] > 0);
1059 //__________________________________________________________________________
1060 void AliITSPreprocessorSPD::PrintChannels(void)
1062 // Print all found bad channels to stdout
1064 // Input : fBadChannelsObjArray
1065 // Output: (dump to stdout)
1070 AliITSChannelSPD *channel = 0;
1072 // Print the bad channels stores in the array
1073 AliInfo("\nModule #\tColumn #\tRow #\n------------------------------------------------");
1074 for (UInt_t module = 0; module < fNumberOfModules; module++)
1077 while (j < fNumberOfBadChannels[module])
1079 channel = (AliITSChannelSPD *) fBadChannelsObjArray->At(i++);
1080 std::cout << module << "\t\t" << channel->GetColumn() << "\t\t" << channel->GetRow() << std::endl;
1082 // Go to next bad channel
1087 AliInfo(Form("%d bad channels were found", fBadChannelsObjArray->GetEntries()));
1091 //__________________________________________________________________________
1092 void AliITSPreprocessorSPD::MarkNoisyChannels(void)
1094 // WARNING: THIS METHOD DOESN'T WORK!!!
1096 // Mark all identified noisy channels
1098 // Input : List of noisy channels, original digits tree
1099 // Output: New digits tree containing SPD digits marked when noisy
1102 // The original digits tree (digitsTree) is cloned except for the SPD branch (ITSDigitSPD).
1103 // This branch is then redefined for each event and will contain all the original
1104 // information. All known noisy channels will be marked by using the TObject status bits
1105 // according to the following scheme. Dead channels are included for completeness. Note
1106 // that a dead channel will NEVER show up among digits..
1108 // Interpretation of digit status bits (LSB):
1109 // Dead channel Noisy channel | Integer
1110 // -----------------------------------------
1115 // meaning e.g. that a channel that is noisy will have the first bit set in its status bits
1117 // Do not continue unless we are processing DAQ data
1120 AliInfo("Marking bad channels");
1122 // Create the storage container that will be used to access the bad channels
1123 if (!fBadChannelsContainer)
1125 // Add the bad channels array to the storage container
1126 // (ownership is passed to the AliRunDataStorage object)
1127 fBadChannelsContainer = new AliITSBadChannelsSPD();
1129 // Convert the bad channels from TObjArray to Int_t[]
1130 AliITSPreprocessorSPD::ConvertObjToIntArray();
1132 // Store the arrays in the bad channels container object
1133 const Int_t kBadChannelsArraySize =
1134 2*fBadChannelsObjArray->GetEntries() + fNumberOfModules;
1135 fBadChannelsContainer->Put(fBadChannelsIntArray, kBadChannelsArraySize,
1136 fBadChannelsIndexArray, fNumberOfModules);
1139 // Create the bad channels helper object
1140 // (will be used to find a bad channel within a TObjArray)
1141 AliITSBadChannelsAuxSPD *aux = new AliITSBadChannelsAuxSPD();
1143 AliITSdigitSPD *digitSPD = 0;
1144 UInt_t numberOfDigits;
1146 Bool_t mark = kFALSE;
1148 TBranch *digitsBranch = 0;
1151 // Create an empty SPD digit array
1152 TObjArray *digitsArraySPD = new TObjArray();
1154 // Get the digits in update mode (we want to modify them if there are noisy channels)
1155 fITSLoader->UnloadDigits();
1156 fITSLoader->LoadDigits("update");
1158 // Get the number of events
1159 UInt_t numberOfEvents = (fRunLoader->TreeE()) ? static_cast<UInt_t>(fRunLoader->TreeE()->GetEntries()) : 0;
1161 // Loop over all events
1162 for (UInt_t event = 0; event < numberOfEvents; event++)
1164 if (event%100 == 0) AliInfo(Form("Event #%d", event));
1166 // Get the current event
1167 fRunLoader->GetEvent(event);
1169 // Get the ITS digits tree
1170 digitsTree = fITSLoader->TreeD();
1172 // Get SPD branch that will contain all digits with marked noisy channels
1173 digitsBranch = digitsTree->GetBranch("ITSDigitsSPD");
1174 digitsBranch->SetAddress(&digitsArraySPD);
1176 // Get the stored number of modules
1177 UInt_t numberOfModules = (Int_t)digitsTree->GetEntries();
1178 TObjArray **newDigitsArraySPD = new TObjArray*[numberOfModules];
1180 Int_t *digitNumber = new Int_t[numberOfModules];
1181 for (UInt_t m = 0; m < numberOfModules; m++)
1183 newDigitsArraySPD[m] = new TObjArray();
1187 AliInfo(Form("ent = %d", (Int_t)digitsTree->GetEntries()));
1189 // Reset the SPD digit arrays to make sure they are empty
1190 digitsArraySPD->Clear();
1192 // Get the SPD digits branch from the original digits tree and set the address
1193 digitsBranch = digitsTree->GetBranch("ITSDigitsSPD");
1194 digitsBranch->SetAddress(&digitsArraySPD);
1196 // Loop over all modules
1197 for (UInt_t module = 0; module < fNumberOfModules; module++)
1199 // Get event data for current module
1200 digitsTree->GetEvent(module);
1202 // Get the hits in the current module
1203 TObjArray *moduleObjArray = fBadChannelsContainer->CreateModuleObjArray(module);
1205 // Get the number of entries
1206 numberOfDigits = digitsArraySPD->GetEntries();
1208 // Loop over all digits and all channels
1209 for (UInt_t digit = 0; digit < numberOfDigits; digit++)
1211 // Get the current digit
1212 digitSPD = (AliITSdigitSPD*) digitsArraySPD->At(digit);
1213 newDigit[0] = digitSPD->GetCoord1(); // column
1214 newDigit[1] = digitSPD->GetCoord2(); // row
1215 newDigit[2] = digitSPD->GetSignal(); // signal
1217 // Check if this channel is noisy
1218 // (Compare with all stored channels in the bad channels array)
1219 if (aux->Find(digitSPD, moduleObjArray))
1221 // Set the mark flag and break the loop
1225 // Store this digit in the SPD digits array using a placement new operation
1226 new ((*newDigitsArraySPD[module])[digitNumber[module]]) AliITSdigitSPD(newDigit);
1228 // Mark it if noisy and store in the noisy channel array
1231 // Store this digit in the marked SPD digits array using a placement new operation
1232 //new ((*badChannels[m])[numberOfBadChannels[m]]) AliITSChannelSPD(newBadChannel);
1233 //new ((*newDigitsArraySPD[module])[digitNumber[module]]) AliITSdigitSPD(newDigit);
1235 // Mark the original channel as noisy
1236 ((*newDigitsArraySPD[module])[digitNumber[module]])->SetBit(kNoisyChannel);
1241 digitNumber[module]++;
1246 delete moduleObjArray;
1249 } // end module loop
1251 digitsBranch->Reset();
1252 digitsBranch->ResetAddress();
1255 delete digitsArraySPD;
1257 digitsTree->Reset();
1259 // WHY THIS RANGE?????????????????????????????????????????????????????????????????????
1260 for (UInt_t n = 0; n < event; n++)
1262 digitsTree->SetBranchAddress("ITSDigitsSPD", &newDigitsArraySPD[n]);
1266 digitsTree->AutoSave();
1269 for (UInt_t n = 0; n < event; n++)
1271 delete newDigitsArraySPD[n];
1273 delete [] newDigitsArraySPD;
1274 newDigitsArraySPD = 0;
1275 delete [] digitNumber;
1280 } // end loop over all events
1282 // Unload the digits
1283 fITSLoader->UnloadDigits();
1293 //__________________________________________________________________________
1294 Bool_t AliITSPreprocessorSPD::Store(AliCDBId &id, AliCDBMetaData *md)
1296 // Store the bad channels object in the calibration database
1297 // (See the corresponding run macro for further explanations)
1299 // Input : fBadChannelsObjArray (now containing all found bad channels), object meta data
1300 // Output: Database file containing the bad channels
1301 // Return: kTRUE if successful
1303 Bool_t status = kFALSE;
1305 AliInfo("Storing bad channels");
1307 // Add the bad channels array to the storage container
1308 // (ownership is passed to the AliRunDataStorage object)
1309 fBadChannelsContainer = new AliITSBadChannelsSPD();
1311 // Convert the bad channels from TObjArray to Int_t[]
1312 AliITSPreprocessorSPD::ConvertObjToIntArray();
1314 // Store the arrays in the bad channels container object
1315 const Int_t kBadChannelsArraySize =
1316 2*fBadChannelsObjArray->GetEntries() + fNumberOfModules;
1317 fBadChannelsContainer->Put(fBadChannelsIntArray, kBadChannelsArraySize,
1318 fBadChannelsIndexArray, fNumberOfModules);
1320 // Store the container
1321 if(!AliCDBManager::Instance()->IsDefaultStorageSet()) {
1322 //AliError("No storage set!");
1324 AliCDBManager::Instance()->SetDefaultStorage("local://Calib");
1327 if (AliCDBManager::Instance()->GetDefaultStorage()->Put(fBadChannelsContainer, id, md))
1329 AliInfo("Bad channels object stored in database");
1334 AliError("Failed to store object in database");
1341 //__________________________________________________________________________
1342 void AliITSPreprocessorSPD::ConvertObjToIntArray()
1344 // Convert the bad channel TObjArray to an Int_t array
1346 // Input : fBadChannelsObjArray (now containing all found bad channels)
1347 // Output: fBadChannelsIntArray
1351 // The TObjArray of this class (fBadChannelsObjArray) is converted to a sequential
1352 // Int_t array (fBadChannelsIntArray) in this method. For each module, the first
1353 // stored number is the number of bad channels in the current module. This is
1354 // followed by all the columns and rows of the bad channels:
1356 // badChannelsArray =
1357 // | N(m) | col0 | row0 | .. | colN(m) | N(m+1) | col0 | row0 | ...
1358 // . .......... module m ......... . .... module m+1 ......
1360 // The bad channels index array (fBadChannelsIndexArray) contains the indices of
1361 // the badChannelsArray, i.e. where the bad channels in certain module starts:
1363 // fBadChannelsObjArray =
1364 // | i0 | i1 | .. | iM | (where M = the number of SPD modules)
1366 // e.g. i1 corresponds to the index of the badChannelsArray where N(1) is stored,
1367 // i.e. the number of bad channels for module 1
1369 const Int_t kBadChannelsArraySize =
1370 2*fBadChannelsObjArray->GetEntries() + fNumberOfModules;
1371 fBadChannelsIntArray = new Int_t[kBadChannelsArraySize]; // Will be deleted in dtor
1372 AliITSChannelSPD *channel = 0;
1377 // Loop over all modules
1378 for (UInt_t module = 0; module < fNumberOfModules; module++)
1380 // Encode the number of bad channels of the current module
1381 fBadChannelsIntArray[k++] = fNumberOfBadChannels[module];
1383 // The columns and rows of the fBadChannelsObjArray will be stored sequentially
1384 // in the Int_t array
1386 while (j < fNumberOfBadChannels[module])
1388 channel = (AliITSChannelSPD *) fBadChannelsObjArray->At(i++);
1389 fBadChannelsIntArray[k++] = channel->GetColumn();
1390 fBadChannelsIntArray[k++] = channel->GetRow();
1392 // Go to next bad channel
1399 //__________________________________________________________________________
1400 Bool_t AliITSPreprocessorSPD::Store(AliCDBId& /*id*/, AliCDBMetaData* /*md*/, Int_t runNumber)
1402 // Store the bad channels object in the calibration database
1403 // (See the corresponding run macro for further explanations)
1405 // Input : fBadChannelsObjArray (now containing all found bad channels), object meta data
1406 // Output: Database file containing the bad channels
1407 // Return: kTRUE if successful
1409 Bool_t status = kFALSE;
1411 AliInfo("Storing bad channels");
1413 AliCDBEntry *entrySPD = AliCDBManager::Instance()->Get("ITS/Calib/CalibSPD", runNumber);
1416 AliWarning("Calibration object retrieval failed! Dummy calibration will be used.");
1418 AliCDBStorage *localStor = AliCDBManager::Instance()->GetStorage("local://$ALICE_ROOT");
1419 entrySPD = localStor->Get("ITS/Calib/CalibSPD", runNumber);
1421 AliFatal("Cannot find SPD calibration entry!");
1426 TObjArray *respSPD = (TObjArray *)entrySPD->GetObject();
1429 AliWarning("Can not get calibration from calibration database !");
1435 AliITSChannelSPD *channel = 0;
1436 AliITSCalibrationSPD* res;
1437 for (Int_t module=0; module<respSPD->GetEntries(); module++) {
1438 res = (AliITSCalibrationSPD*) respSPD->At(module);
1440 while (j < fNumberOfBadChannels[module]) {
1441 channel = (AliITSChannelSPD *) fBadChannelsObjArray->At(i++);
1442 res->AddDead(channel->GetColumn(),channel->GetRow());
1443 // Go to next bad channel
1449 AliCDBStorage *storage = AliCDBManager::Instance()->GetDefaultStorage();
1451 AliWarning("No default storage set! Will use dummy one");
1452 storage = AliCDBManager::Instance()->GetStorage("local://$ALICE_ROOT");
1453 if(!storage) AliFatal("Could not even set dummy storage! Something very strange is happening...");
1456 status = storage->Put(entrySPD);
1457 entrySPD->SetObject(NULL);
1458 entrySPD->SetOwner(kTRUE);