Calibration framework improved (E. Crescio)

[u/mrichter/AliRoot.git] / ITS / AliITSPreprocessorSPD.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/*
$Log$
Revision 1.1  2005/10/11 12:31:50  masera
Preprocessor classes for SPD (Paul Nilsson)

*/

///////////////////////////////////////////////////////////////////////////
// AliITSPreprocessorSPD implementation by P. Nilsson 2005
// MAIN AUTHOR/CONTACT: Paul.Nilsson@cern.ch
//
// The purpose of this class is to provide algorithms for identification
// of "bad channels" such as dead channels and noisy channels in the SPD
//
// Examples on how to use this class can be found in the macros:
//
// .findNoisyChannels.C (Locate and store noisy channels in the CDB)
// .readNoisyChannels.C (Read noisy channels from the CDB)
// .findDeadChannels.C  (Locate and store dead channels in the CDB)
// .readDeadChannels.C  (Read dead channels from the CDB)
///////////////////////////////////////////////////////////////////////////

#include "TFile.h"
#include "AliITSPreprocessorSPD.h"
#include "AliCDBEntry.h"
#include "AliITSCalibrationSPD.h" 
ClassImp(AliITSPreprocessorSPD)


//__________________________________________________________________________
AliITSPreprocessorSPD::AliITSPreprocessorSPD(void):
  fITSLoader(0x0),
  fRunLoader(0x0),
  fThresholdRatio(5.),
  fThreshold(5),
  fMaximumNumberOfEvents(1000000),
  fHighestModuleNumber(0),
  fSelectedAlgorithm(kOptimizedForRealData),
  fGeometryMode(kALICEGeometry),
  fNumberOfBadChannels(0),
  fInit(kFALSE),
  fVMEMode(kFALSE),
  fDigitsHistogram(0),
  fBadChannelsObjArray(0),
  fBadChannelsIntArray(0),
  fBadChannelsIndexArray(0),
  fBadChannelsContainer(0)
{
  // Default constructor for the SPD preprocessor
  //
  // Initialization has to be done by hand using Set* methods and the Open method
  //
  // Input :
  // Output:
  // Return: <empty/uninitialized AliITSPreprocessorSPD object>
}


//__________________________________________________________________________
AliITSPreprocessorSPD::AliITSPreprocessorSPD(const char *fileName, const char *mode,
                                             const char *fileNameg, const Int_t maxNumberOfEvents):
  fITSLoader(0x0),
  fRunLoader(0x0),
  fInit(kFALSE)
{
  // Standard constructor for the SPD preprocessor
  //
  // Input : Name of digit file
  // Output: fInit
  // Return: Initialized SPD preprocessor object

  // Initialize
  AliITSPreprocessorSPD::Init();
  AliITSPreprocessorSPD::SetMaximumNumberOfEvents(maxNumberOfEvents);

  // Open and read the galice and digit files
  if (!AliITSPreprocessorSPD::Open(fileName, mode, fileNameg))
    {
      AliError("Failed to open file");
      fInit = kFALSE;
    }
}


//__________________________________________________________________________
AliITSPreprocessorSPD::AliITSPreprocessorSPD(const AliITSPreprocessorSPD &prep) :
  TTask(prep)
{
  // Default copy constructor
  // Notice that only pointer addresses are copied!
  // Memory allocations of new objects are not done.

  fITSLoader = prep.fITSLoader;
  fRunLoader = prep.fRunLoader;
  fThresholdRatio = prep.fThresholdRatio;
  fThreshold = prep.fThreshold;
  fMaximumNumberOfEvents = prep.fMaximumNumberOfEvents;
  fHighestModuleNumber = prep.fHighestModuleNumber;
  fSelectedAlgorithm = prep.fSelectedAlgorithm;
  fGeometryMode = prep.fGeometryMode;
  fNumberOfBadChannels = prep.fNumberOfBadChannels;
  fInit = prep.fInit;
  fVMEMode = prep.fVMEMode;
  fDigitsHistogram = prep.fDigitsHistogram;
  fBadChannelsObjArray = prep.fBadChannelsObjArray;
  fBadChannelsIntArray = prep.fBadChannelsIntArray;
  fBadChannelsIndexArray = prep.fBadChannelsIndexArray;
  fBadChannelsContainer = prep.fBadChannelsContainer;
}


//__________________________________________________________________________
AliITSPreprocessorSPD& AliITSPreprocessorSPD::operator=(const AliITSPreprocessorSPD &prep)
{
  // Assignment operator

  AliError("Not implemented");

  if (this != &prep) { } // Do not delete this line

  return *this;
}


//__________________________________________________________________________
AliITSPreprocessorSPD::~AliITSPreprocessorSPD(void)
{
  // Destructor for the SPD preprocessor

  if (fDigitsHistogram)
    {
      // try fDigitsHistogram->Delete(); if the following crashes
      for (UInt_t module = 0; module < fNumberOfModules; module++)
        {
          (*fDigitsHistogram)[module]->Delete();
        }
      delete fDigitsHistogram;
      fDigitsHistogram = 0;
    }

  if (fNumberOfBadChannels)
    {
      delete [] fNumberOfBadChannels;
      fNumberOfBadChannels = 0;
    }

  if (fBadChannelsIntArray)
    {
      delete [] fBadChannelsIntArray;
      fBadChannelsIntArray = 0;
    }

  if (fBadChannelsIndexArray)
    {
      delete [] fBadChannelsIndexArray;
      fBadChannelsIndexArray = 0;
    }

  if (fBadChannelsObjArray)
    {
      fBadChannelsObjArray->Delete();
      delete fBadChannelsObjArray;
      fBadChannelsObjArray = 0;
    }

  delete fRunLoader;
  fRunLoader = 0;

  delete fBadChannelsContainer;
  fBadChannelsContainer = 0;
}


//__________________________________________________________________________
void AliITSPreprocessorSPD::Init(void)
{
  // Initialization of the SPD preprocessor
  //
  // Input : (void)
  // Output: Several logistical variables
  // Return: (void)

  // Default maximum number of events per histogram
  fMaximumNumberOfEvents = 1000000;

  // Default noisy channel removal algorithm
  fSelectedAlgorithm = kOptimizedForRealData;

  // Default noisy channel threshold and threshold ratio
  // (threshold for current bin content divided by the average neighboring bin contents)
  fThreshold = 10;
  fThresholdRatio = 5.;

  // Set default geometry mode (ALICE geometry). This also sets fNumberOfModules
  AliITSPreprocessorSPD::SetGeometryMode(kALICEGeometry);

  // The highest module number with found bad channels
  fHighestModuleNumber = 0;

  // Assume input is not from a VME file
  fVMEMode = kFALSE;

  // Initialization is complete
  fInit = kTRUE;
}


//__________________________________________________________________________
void AliITSPreprocessorSPD::SetGeometryMode(UInt_t mode)
{
  // Set the geometry mode
  //
  // Input : Geometry mode (either kTestBeamGeometry or kALICEGeometry)
  // Output: 
  // Return: (void)

  fGeometryMode = mode;

  // In case of an input VME file, the number of modules has already been fixed.
  // Do not try to change it
  if (!fVMEMode)
    {
      if (mode == kALICEGeometry)
        {
          fNumberOfModules = kNumberOfSPDModules;
        }
      else if (mode == kTestBeamGeometry)
        {
          fNumberOfModules = kNumberOfTestBeamSPDModules;
        }
      else
        {
          AliError("Unknown geometry mode, defaults to ALICE geometry");
          fNumberOfModules = kNumberOfSPDModules;
        }
    }
}


//__________________________________________________________________________
void AliITSPreprocessorSPD::SetFakeNoisyChannel(Int_t module, Int_t column, Int_t row)
{
  // Introduce a fake noisy channel in the hit histograms
  //
  // Input : Module, column and row numbers
  // Output: Updated hit histograms
  // Return: (void)

  if ((UInt_t)module < fNumberOfModules)
    {
      ((TH2F*)(*fDigitsHistogram)[module])->Fill(column, row, 1000000);
    }
  else
    {
      AliError("No such module number");
    }
}


//__________________________________________________________________________
Bool_t AliITSPreprocessorSPD::Open(const char *fileName, const char *mode, const char *fileNameg)
{
  // Open digit file
  //
  // Input : Name and mode of ITS digit file, name of galice file
  // Output: Digits
  // Return: kTRUE if loaders succeed

  Bool_t status = kFALSE;
  Bool_t status0 = kFALSE;
  Bool_t status1 = kFALSE;
  Bool_t status2 = kFALSE;
  Bool_t status3 = kFALSE;

  // Only proceed if initialization has been done
  if (fInit)
    {
      TString m = (TString)mode;
      if (m == "daq" || m == "DAQ")
        {
          // Open the data file and get the run loader
          fRunLoader = AliRunLoader::Open(fileNameg);
          if (fRunLoader)
            {
              // Get the gAlice object
              status0 = AliITSPreprocessorSPD::GetgAlice();

              // Get the ITS digits
              if (status0) status1 = AliITSPreprocessorSPD::GetITSDigits(fileName);

              // Create the test beam object
              if (status1) status2 = AliITSPreprocessorSPD::CreateGeometryObj();

              // Fill histograms with DAQ digits
              if (status2) status3 = AliITSPreprocessorSPD::FillHistograms();

              status = status0 & status1 & status2 & status3;
            }
          else
            {
              AliError("Failed to get the run loader");
            }
        }
      else if (m == "vme" || m == "VME")
        {
          // Open the VME file that contains histograms with fired channels as read by the stand-alone VME system
          TFile *vmeFile = TFile::Open(fileName);

          if (!vmeFile->IsOpen())
            {
              AliError("Could not open VME input file");
            }
          else
            {
              // Get the histograms from the VME file that contains all fired channels
              status0 = AliITSPreprocessorSPD::GetVMEHistograms(vmeFile);

              // Create the test beam object
              if (status0) status1 = AliITSPreprocessorSPD::CreateGeometryObj();
            }

          // This boolean will be used to override any attempts of changing the number of modules by the user
          // with the SetGeometryMode method. For VME files the number of modules will entirely be determined by
          // the input VME root file, i.e. by the number of histograms in this file
          fVMEMode = kTRUE;

          status = status0 & status1;
        }
      else
        {
          AliError("Unknown filetype - assuming DAQ file");
        }

      // At this stage, the final number of modules will be known (safe to define arrays)
      // In case data is read from a VME root file, it will not be known before
      if (status)
        {
          // Create the arrays for bookkeeping and storing the noisy channels
          if (!fBadChannelsObjArray)
            {
              fBadChannelsObjArray = new TObjArray();
            }
          if (!fBadChannelsIndexArray)
            {
              // This array will contain the begin and end indices for each module, i.e. where to begin
              // and stop reading the fBadChannelsObjArray for a certain module.
              // The last position of the array will contain the end index of the last module
              fBadChannelsIndexArray = new Int_t[fNumberOfModules + 1];
              for (UInt_t module = 0; module < fNumberOfModules + 1; module++)
                {
                  fBadChannelsIndexArray[module] = 0;
                }
            }
        }
    }
  else
    {
      AliError("SPD preprocessor not initialized. Can't load digits");
    }

  return status;
}


//__________________________________________________________________________
Bool_t AliITSPreprocessorSPD::GetgAlice(void)
{
  // Get the gAlice object
  //
  // Input : (void)
  // Output: gAlice object
  // Return: kTRUE if the gAlice object was found

  Bool_t status = kTRUE;

  // Load gAlice
  fRunLoader->LoadgAlice();
  if (!fRunLoader->GetAliRun())
    {
      AliError("gAlice not found on file. Aborting.");
      status = kFALSE;
    }

  return status;
}


//__________________________________________________________________________
Bool_t AliITSPreprocessorSPD::GetVMEHistograms(TFile *vmeFile)
{
  // Get pre-filled digit histograms from input VME file
  //
  // Input : pointer to VME file
  // Output: 
  // Return: kTRUE if histograms are found on file

  Bool_t status = kFALSE;

  // Get the file directory
  TDirectory *dir = (TDirectory *)vmeFile;

  // Get the number of keys (histograms in this case corresponding to modules/ladders)
  fNumberOfModules = dir->GetNkeys();
  if ((fNumberOfModules > 0) && (fNumberOfModules <= kNumberOfSPDModules))
    {
      status = kTRUE;

      // Create bad channel histograms
      fDigitsHistogram = new TObjArray(fNumberOfModules);

      // Create a key iterator
      TIter nextkey(dir->GetListOfKeys());
      TKey *key = 0;

      // Loop over all objects, read them in to memory one by one
      UInt_t module = 0;
      while ( (key = (TKey *)nextkey()) )
        {
          (*fDigitsHistogram)[module++] = (TH2F *)key->ReadObj();

          // For safety
          if (module > fNumberOfModules)
            {
              AliError("Wrong number of keys in VME file");
              status = kFALSE;
              break;
            }
        }
    }
  else
    {
      AliError("Wrong number of histograms in VME file");
    }
  return status;
}

//__________________________________________________________________________
Bool_t AliITSPreprocessorSPD::GetITSDigits(const char *fileName)
{
  // Get the ITS digits
  //
  // Input : name of digits file
  // Output: fITSLoader object, ITS digits
  // Return: kTRUE if loader succeed

  Bool_t status = kTRUE;

  // Load the gAlice and the header
  fRunLoader->LoadgAlice();
  fRunLoader->LoadHeader();

  // Get the ITS loader
  fITSLoader = (AliITSLoader*) fRunLoader->GetLoader("ITSLoader");
  if (!fITSLoader)
    {
      AliError("ITS loader not found");
      status = kFALSE;
    }
  else
    {
      // Open the digits file
      fITSLoader->SetDigitsFileName(fileName);
    }

  return status;
}


//__________________________________________________________________________
TClonesArray* AliITSPreprocessorSPD::CreateDigitsArray(void) const
{
  // Creation of the digits array
  //
  // Input : (void)
  // Output:
  // Return: Pointer to the SPD digits array

  // Create an array for 5 chips with 8192 channels each
  TClonesArray *digitsArray = new TClonesArray("AliITSdigitSPD", kNumberOfChannels);

  return digitsArray;
}


//__________________________________________________________________________
Bool_t AliITSPreprocessorSPD::CreateGeometryObj(void)
{
  // Creation of the geometry object
  //
  // This object is used to get the number of SPD half-staves
  //
  // Input : (void)
  // Output: fGeometryObj
  // Return: kTRUE if fGeometryObj has been created

  Bool_t status = true;

  // Create geometry object
  // fGeometryObj = new ...
  // if (!fGeometryObj) status = kFALSE;

  // Get SPD parameters
  // fNumberOfColumns = fGeometryObject->GetNumberOfColumns();
  // fNumberOfRows = fGeometryObject->GetNumberOfRows();

  fNumberOfColumns = kNumberOfColumns;
  fNumberOfRows = kNumberOfRows;

  return status;
}


//__________________________________________________________________________
void AliITSPreprocessorSPD::CreateHistograms(void)
{
  // Create the noisy channel histograms
  //
  // Input : (void)
  // Output: Noisy channel histograms
  // Return: (void)

  TString moduleName;
  char n[4]; // For storing the module number (maximum module number is 240)

  // Create noisy channel histograms
  fDigitsHistogram = new TObjArray(fNumberOfModules);

  for (UInt_t i = 0; i < fNumberOfModules; i++)
    {
      // Histogram name
      moduleName = "module_";
      sprintf(n,"%d",i);
      moduleName.Append(n);

      (*fDigitsHistogram)[i] = new TH2F(moduleName,"Digits",
                                        fNumberOfColumns,-.5,(1.*fNumberOfColumns - .5),
                                        fNumberOfRows,-.5,(1.*fNumberOfRows - .5));
    }
}


//__________________________________________________________________________
Bool_t AliITSPreprocessorSPD::FillHistograms(void)
{
  // Fill the histograms with digits (hit positions of unclustered hits)
  //
  // (There is one digit histogram per SPD module, i.e. a half-stave, 10 chips)
  //
  // Input : No arguments (Empty digit histograms)
  // Output: Filled digit histograms
  // Return: kTRUE if histograms are filled with digits, kFALSE otherwise

  AliInfo("Filling noisy channel histograms");

  Bool_t status = kTRUE;
  AliITSdigitSPD *digitSPD = 0;
  UInt_t row;
  UInt_t column;
  TBranch *digitsSPDBranch;
  TTree *digitsTree;

  // Get the digits
  fITSLoader->LoadDigits("read");

  // Create noisy channel histograms
  AliITSPreprocessorSPD::CreateHistograms();

  // Create an empty SPD digits array
  TClonesArray *digitsArraySPD = AliITSPreprocessorSPD::CreateDigitsArray();

  // Get number of events
  UInt_t numberOfEvents = (fRunLoader->TreeE()) ? static_cast<UInt_t>(fRunLoader->TreeE()->GetEntries()) : 0;

  // Make sure we don't try to analyze more data than there actually is
  if (numberOfEvents < fMaximumNumberOfEvents)
    {
      fMaximumNumberOfEvents = numberOfEvents;
    }

  // Loop over all digits and put them in the corresponding histograms
  if (numberOfEvents > 0)
    {
      for (UInt_t event = 0; event < fMaximumNumberOfEvents; event++)
        {
          // Get the current event
          fRunLoader->GetEvent(event);

          // Get the ITS digits tree
          digitsTree = fITSLoader->TreeD();

          // Disable all branches except the SPD branch to speed up the reading process
          digitsTree->SetBranchStatus("ITSDigitsSPD",1);
          digitsTree->SetBranchStatus("ITSDigitsSDD",0);
          digitsTree->SetBranchStatus("ITSDigitsSSD",0);

          // Reset the SPD digits array to make sure it is empty
          digitsArraySPD->Clear();

          // Get the SPD digits branch and set the address
          digitsSPDBranch = digitsTree->GetBranch("ITSDigitsSPD");

          digitsSPDBranch->SetAddress(&digitsArraySPD);

          if (event%100 == 0) AliInfo(Form("Event #%d", event));

          // Loop over all modules
          UInt_t numberOfDigits = 0;
          for (UInt_t module = 0; module < fNumberOfModules; module++)
            {
              // Get event data for current module
              digitsTree->GetEvent(module);

              // Get the number of entries
              numberOfDigits = digitsArraySPD->GetEntries();

              // Loop over all digits
              for (UInt_t digit = 0; digit < numberOfDigits; digit++)
                {
                  // Get the current digit
                  digitSPD = (AliITSdigitSPD*) digitsArraySPD->At(digit);
                  row = digitSPD->GetCoord1();
                  column = digitSPD->GetCoord2();

                  // Fill the digits histogram
                  ((TH2F*)(*fDigitsHistogram)[module])->Fill(column, row);
                  
                } // end digit loop
            } // end module loop
        } // end event loop
    } // end numberOfEvents > 0

  else
    {
      status = kFALSE;
    }

  // Cleanup
  delete digitsArraySPD;
  digitsArraySPD = 0;

  return status;
}


//__________________________________________________________________________
Bool_t AliITSPreprocessorSPD::FindDeadChannels(void)
{
  // Locate dead channels
  //
  // Input : Filled hit histograms
  // Output: TObjArray (fBadChannelsObjArray) with all identified dead channels
  // Return: kTRUE if dead channels have been found

  Bool_t status = kFALSE;

  // Proceed only if properly initialized
  if (fInit)
    {
      if (fVMEMode)
        {
          // Initialize counters
          fNumberOfBadChannels = new Int_t[fNumberOfModules];
          for (UInt_t module = 0; module < fNumberOfModules; module++)
            {
              fNumberOfBadChannels[module] = 0;
            }

          // Loop over all modules (intentional modularization - DCS will occationally want to
          // look for noisy channels in certain modules, but not all)
          fIndex = 0; // Global bad channels array counter (must be reset here)
          for (UInt_t module = 0; module < fNumberOfModules; module++)
            {
              status |= AliITSPreprocessorSPD::FindDeadChannelsInModule(module);
            }
        }
      else
        {
          AliError("Dead channels can only be found in data taken with stand-alone VME system");
        }
    }
  else
    {
      AliError("Not properly initialized");
    }

  return status;
}


//__________________________________________________________________________
Bool_t AliITSPreprocessorSPD::FindDeadChannelsInModule(UInt_t module)
{
  // Locate dead channels
  // This method assumes that the preprocessor is operator in VME mode.
  // The algorithm is very simple. It assumes that the data was taken in
  // a mode where all working SPD pixels should respond as being hit.
  // fThreshold is used as the limit where everything below this value will
  // be considered as "dead".
  //
  // Input : Filled hit histograms
  // Output: TObjArray (fBadChannelsObjArray) with all identified noisy channels
  // Return: kTRUE if dead channels have been found

  // Store the index number for this module
  fBadChannelsIndexArray[module] = fIndex++;

  UInt_t xBin, numberOfXBins;
  UInt_t yBin, numberOfYBins;
  Float_t binContent;

  numberOfXBins = ((TH2F*)(*fDigitsHistogram)[module])->GetNbinsX();
  numberOfYBins = ((TH2F*)(*fDigitsHistogram)[module])->GetNbinsY();

  // Loop over all bins in this histogram
  for (xBin = 1; xBin <= numberOfXBins; xBin++)
    for (yBin = 1; yBin <= numberOfYBins; yBin++)
      {
        binContent = ((TH2F*)(*fDigitsHistogram)[module])->GetBinContent(xBin, yBin);

        // Store this channel/bin if outside accepted limits
        // A channel has to fire MORE times than the fThreshold value, or it will
        // be considered as "dead"
        if (binContent <= fThreshold)
          {
            // Store the dead channel in the array
            // The channel object will be deleted in the destructor using the TObjArray Delete() method
            // (The array will assume ownership of the channel object)
            AliITSChannelSPD *channel = new AliITSChannelSPD((Int_t)(xBin - 1), (Int_t)(yBin - 1));

            // Store the noisy channel in the array
            fBadChannelsObjArray->Add(channel);

            // Keep track of the number of bad channels in this module
            fNumberOfBadChannels[module]++;
            fIndex += 2;

            // Keep track of the highest module number
            if (module > fHighestModuleNumber) fHighestModuleNumber = module;

            //AliInfo(Form("New dead pixel in (m,c,r) = (%d,%d,%d)", module, xBin - 1, yBin - 1));
          }
      } // end bin loop

  return (fNumberOfBadChannels[module] > 0);
}


//__________________________________________________________________________
Bool_t AliITSPreprocessorSPD::FindNoisyChannels(void)
{
  // Locate noisy channels by searching through the digit histograms
  // (There is one digit histogram per SPD module, i.e. a half-stave, 5 chips)
  //
  // Input : Digits
  // Output: TObjArray (fBadChannelsObjArray) with all identified noisy channels
  // Return: kTRUE if noisy channels have been found

  Bool_t status = kFALSE;

  // Proceed only if properly initialized
  if (fInit)
    {
      // (For testing purposes, noisy channels can be inserted here)
      //SetFakeNoisyChannel(4,10,10);
      // Initialize counters
      fNumberOfBadChannels = new Int_t[fNumberOfModules];
      for (UInt_t module = 0; module < fNumberOfModules; module++)
        {
          fNumberOfBadChannels[module] = 0;
        }

      // Scan through all the histogram bins and search for average filling deviations
      if ((fSelectedAlgorithm == kOptimizedForRealData) || (fSelectedAlgorithm == kOptimizedForRealDataRMS))
        {
          // Noisy channel algorithm optimized for real data..........................
          // Histograms can have any shape (both thresholds and quotients are used)
          // This algorithm can be used to find noisy channels even if the data was
          // taken with beam. All channels outside accepted limits (set by fThreshold
          // and fThresholdRatio) will be marked as noisy

          if (fSelectedAlgorithm == kOptimizedForRealData)
            {
              AliInfo("Searching for noisy channels (optimized for real data)");
            }
          else
            {
              AliInfo("Searching for noisy channels (optimized for real data, RMS version)");
            }

          // Loop over all modules (intentional modularization - DCS will occationally want to
          // look for noisy channels in certain modules, but not all)
          fIndex = 0; // Global bad channels array counter (must be reset here)
          for (UInt_t module = 0; module < fNumberOfModules; module++)
            {
              status |= AliITSPreprocessorSPD::FindNoisyChannelsInModuleAlgo0(module);
            }
        } // end algorithm 0
      else
        {
          // Noisy channel algorithm optimized for calibration data...........................
          // Histograms will/should only contain noisy channels (only thresholds are used)
          // Calibration histograms should have no background. The calibration data
          // should have been taken without beam. All channels outside accepted limits
          // (set by fThreshold) will be marked as noisy

          AliInfo("Searching for noisy channels (optimized for calibration data)");

          // Loop over all modules (intentional modularization - DCS will occationally want to
          // look for noisy channels in certain modules, but not all)
          fIndex = 0; // Global bad channels array counter (must be reset here)
          for (UInt_t module = 0; module < fNumberOfModules; module++)
            {
              status |= AliITSPreprocessorSPD::FindNoisyChannelsInModuleAlgo1(module);
            }
        } // end algorithm 1
    } // end if fInit
  else
    {
      AliError("Not properly initialized");
    }

  return status;
}


//__________________________________________________________________________
Bool_t AliITSPreprocessorSPD::FindNoisyChannelsInModuleAlgo0(UInt_t module)
{
  // Locate the noisy channels in a module (optimized for real data)
  //
  // For each channel in these histograms, the algorithm checks the average
  // filling of the neighboring 3, 5 or 8 channels (depending on the location
  // of the channel in question; corner, border or inside), or compares with the
  // RMS of the neighbors. If the average is "much" lower, the channel will be
  // considered to be noisy. The default noisy-to-normal fraction is stored in the
  // fThresholdRatio varible. It can be set with the SetThresholdRatio method.
  // The channel also has to have fired more than a certain minimum, fThreshold.
  // It can be set with the SetThreshold method.
  //
  // To avoid difficulties with noisy channels that occur in pairs, the
  // neighboring channel with largest number of fillings will be removed from
  // the average calculation.
  //
  // NOTE: Since this method modifies the fBadChannelsObjArray and fBadChannelsIndexArray
  //       it is essential to initialize the fIndex counter before calling this module
  //       the first time. The bad channel data does not have to be ordered per module
  //       in the fBadChannelsObjArray, but the indices of where the data of a certain module
  //       starts has to be correct. A wrong fIndex can lead to segmentation violation
  //
  // Input : module number, filled digit histograms
  // Output: TObjArray (fBadChannelsObjArray) with all identified noisy channels,
  //         Int_t[] (fBadChannelsIndexArray) with fBadChannelsObjArray module indices,
  //         number of noisy channels in this module (global variable fNumberOfBadChannels[module])
  // Return: kTRUE if there are noisy channels in this module

  // Store the index number for this module
  fBadChannelsIndexArray[module] = fIndex++;

  UInt_t xBin, numberOfXBins;
  UInt_t yBin, numberOfYBins;
  UInt_t neighborXBin;
  UInt_t neighborYBin;
  UInt_t numberOfNeighboringBins;
  Float_t binContent;
  Float_t sumBinContent;
  Float_t neighborBinContent;
  Float_t maxBinContent;
  Float_t averageBinContent;
  Float_t ratio;

  numberOfXBins = (UInt_t)((TH2F*)(*fDigitsHistogram)[module])->GetNbinsX();
  numberOfYBins = (UInt_t)((TH2F*)(*fDigitsHistogram)[module])->GetNbinsY();

  // Loop over all bins in this histogram
  for (xBin = 1; xBin <= numberOfXBins; xBin++)
    for (yBin = 1; yBin <= numberOfYBins; yBin++)
      {
        numberOfNeighboringBins = 0;
        averageBinContent = 0.;
        sumBinContent = 0.;
        binContent = ((TH2F*)(*fDigitsHistogram)[module])->GetBinContent(xBin, yBin);
        maxBinContent = 0.;

        // Calculate the average pixel level on the surrounding pixels
        for (neighborXBin = xBin - 1; neighborXBin <= xBin + 1; neighborXBin++)
          for (neighborYBin = yBin - 1; neighborYBin <= yBin + 1; neighborYBin++)
            {
              if ( !((neighborXBin == xBin) && (neighborYBin == yBin)) )
                {
                  // Only update the number of neighboring bins when we are not on the border
                  if ((neighborXBin > 0) && (neighborXBin <= numberOfXBins+1) &&
                      (neighborYBin > 0) && (neighborYBin <= numberOfYBins+1))
                    {
                      neighborBinContent =
                        ((TH2F*)(*fDigitsHistogram)[module])->GetBinContent(neighborXBin, neighborYBin);

                      if (fSelectedAlgorithm == kOptimizedForRealDataRMS)
                        {
                          // RMS
                          sumBinContent += neighborBinContent*neighborBinContent;
                        }
                      else
                        {
                          // Geometrical mean
                          sumBinContent += neighborBinContent;
                        }

                      if (neighborBinContent > maxBinContent) maxBinContent = neighborBinContent;

                      numberOfNeighboringBins++;
                    }
                }
            }
        
        // Calculate the average. Remove the largest neighboring bin
        // (Correction for potential clusters of noisy channels)
        if (fSelectedAlgorithm == kOptimizedForRealDataRMS)
          {
            // Square the max bin content before removing it from the average calculation
            maxBinContent *= maxBinContent;
            // RMS
            averageBinContent = TMath::Sqrt((sumBinContent - maxBinContent)/(Float_t)(numberOfNeighboringBins - 1));
          }
        else
          {
            // Geometrical mean
            averageBinContent = (sumBinContent - maxBinContent)/(Float_t)(numberOfNeighboringBins - 1);
          }

        // Store this channel/bin if outside accepted limits
        // The threshold ratio is the threshold for the current bin content divided by the
        // average neighboring bin contents. The threshold bin content is the minimum number of
        // times a channel has to have fired to be called noisy
        ratio = (averageBinContent > 0) ? binContent/averageBinContent : 0.;
        if ( ((ratio >= fThresholdRatio) || (ratio == 0.)) && (binContent >= fThreshold) )
          {
            // Store the noisy channel in the array
            // The channel object will be deleted in the destructor using the TObjArray Delete() method
            // (The array will assume ownership of the channel object)
            AliITSChannelSPD *channel = new AliITSChannelSPD((Int_t)(xBin - 1), (Int_t)(yBin - 1));

            // Store the noisy channel in the array
            fBadChannelsObjArray->Add(channel);

            // Keep track of the number of bad channels in this module
            fNumberOfBadChannels[module]++;
            fIndex += 2;

            // Keep track of the highest module number
            if (module > fHighestModuleNumber) fHighestModuleNumber = module;

            AliInfo(Form("New noisy pixel in (m,c,r) = (%d,%d,%d)", module, xBin - 1, yBin - 1));
            AliInfo(Form("- Noisy pixel fired %d times, average neighborhood: %f",(Int_t)binContent,averageBinContent));
          }
      } // end bin loop

  return (fNumberOfBadChannels[module] > 0);
}


//__________________________________________________________________________
Bool_t AliITSPreprocessorSPD::FindNoisyChannelsInModuleAlgo1(UInt_t module)
{
  // Locate the noisy channels in a module (optimized for calibration data)
  //
  // This algorithm locates noisy channels by assuming original data was taken
  // in calibration mode. This should be done without beam and will thus only
  // contain data corresponding to background and noisy channels. The latter
  // should be clearly visible in this data so this algorithm simply assumes
  // that all histogram bins that are filled more than fThreshold times are
  // noisy.
  //
  // NOTE: Since this method modifies the fBadChannelsObjArray and fBadChannelsIndexArray
  //       it is essential to initialize the fIndex counter before calling this module
  //       the first time. The bad channel data does not have to be ordered per module
  //       in the fBadChannelsObjArray, but the indices of where the data of a certain module
  //       starts has to be correct. A wrong fIndex can lead to segmentation violation
  //
  // Input : module number, filled digit histograms
  // Output: TObjArray (fBadChannelsObjArray) with all identified noisy channels,
  //         Int_t[] (fBadChannelsIndexArray) with fBadChannelsObjArray module indices,
  //         number of noisy channels in this module (global variable fNumberOfBadChannels[module])
  // Return: kTRUE if there are noisy channels in this module

  // Store the index number for this module
  fBadChannelsIndexArray[module] = fIndex++;

  UInt_t xBin, numberOfXBins;
  UInt_t yBin, numberOfYBins;
  Float_t binContent;

  numberOfXBins = ((TH2F*)(*fDigitsHistogram)[module])->GetNbinsX();
  numberOfYBins = ((TH2F*)(*fDigitsHistogram)[module])->GetNbinsY();

  // Loop over all bins in this histogram
  for (xBin = 1; xBin <= numberOfXBins; xBin++)
    for (yBin = 1; yBin <= numberOfYBins; yBin++)
      {
        binContent = ((TH2F*)(*fDigitsHistogram)[module])->GetBinContent(xBin, yBin);

        // Store this channel/bin if outside accepted limits
        // The threshold bin content is the minimum number of times a channel has to have
        // fired to be called noisy
        if (binContent >= fThreshold)
          {
            // Store the noisy channel in the array
            // The channel object will be deleted in the destructor using the TObjArray Delete() method
            // (The array will assume ownership of the channel object)
            AliITSChannelSPD *channel = new AliITSChannelSPD((Int_t)(xBin - 1), (Int_t)(yBin - 1));

            // Store the noisy channel in the array
            fBadChannelsObjArray->Add(channel);

            // Keep track of the number of bad channels in this module
            fNumberOfBadChannels[module]++;
            fIndex += 2;

            // Keep track of the highest module number
            if (module > fHighestModuleNumber) fHighestModuleNumber = module;

            AliInfo(Form("New noisy pixel in (m,c,r) = (%d,%d,%d)", module, xBin - 1, yBin - 1));
            AliInfo(Form("- Noisy pixel fired %d times",(Int_t)binContent));
          }
      } // end bin loop

  return (fNumberOfBadChannels[module] > 0);
}


//__________________________________________________________________________
void AliITSPreprocessorSPD::PrintChannels(void)
{
  // Print all found bad channels to stdout
  //
  // Input : fBadChannelsObjArray
  // Output: (dump to stdout)
  // Return: (void)

  Int_t i = 0;
  Int_t j = 0;
  AliITSChannelSPD *channel = 0;

  // Print the bad channels stores in the array
  AliInfo("\nModule #\tColumn #\tRow #\n------------------------------------------------");
  for (UInt_t module = 0; module < fNumberOfModules; module++)
    {
      j = 0;
      while (j < fNumberOfBadChannels[module])
        {
          channel = (AliITSChannelSPD *) fBadChannelsObjArray->At(i++);
          std::cout << module << "\t\t" << channel->GetColumn() << "\t\t" << channel->GetRow() << std::endl;

          // Go to next bad channel
          j++;
        }
    }

  AliInfo(Form("%d bad channels were found", fBadChannelsObjArray->GetEntries()));
}


//__________________________________________________________________________
void AliITSPreprocessorSPD::MarkNoisyChannels(void)
{
  // WARNING: THIS METHOD DOESN'T WORK!!!
  //
  // Mark all identified noisy channels
  //
  // Input : List of noisy channels, original digits tree
  // Output: New digits tree containing SPD digits marked when noisy
  // Return: (void)
  //
  // The original digits tree (digitsTree) is cloned except for the SPD branch (ITSDigitSPD).
  // This branch is then redefined for each event and will contain all the original
  // information. All known noisy channels will be marked by using the TObject status bits
  // according to the following scheme. Dead channels are included for completeness. Note
  // that a dead channel will NEVER show up among digits..
  //
  // Interpretation of digit status bits (LSB):
  // Dead channel   Noisy channel   |  Integer
  // -----------------------------------------
  //      0               0         |     0
  //      0               1         |     1
  //      1               0         |     2
  //
  // meaning e.g. that a channel that is noisy will have the first bit set in its status bits

  // Do not continue unless we are processing DAQ data
  if (!fVMEMode)
    {
      AliInfo("Marking bad channels");

      // Create the storage container that will be used to access the bad channels
      if (!fBadChannelsContainer)
        {
          // Add the bad channels array to the storage container
          // (ownership is passed to the AliRunDataStorage object)
          fBadChannelsContainer = new AliITSBadChannelsSPD();

          // Convert the bad channels from TObjArray to Int_t[]
          AliITSPreprocessorSPD::ConvertObjToIntArray();

          // Store the arrays in the bad channels container object
          const Int_t kBadChannelsArraySize =
            2*fBadChannelsObjArray->GetEntries() + fNumberOfModules;
          fBadChannelsContainer->Put(fBadChannelsIntArray, kBadChannelsArraySize,
                                     fBadChannelsIndexArray, fNumberOfModules);
        }

      // Create the bad channels helper object
      // (will be used to find a bad channel within a TObjArray)
      AliITSBadChannelsAuxSPD *aux = new AliITSBadChannelsAuxSPD();

      AliITSdigitSPD *digitSPD = 0;
      UInt_t numberOfDigits;
      Int_t newDigit[3];
      Bool_t mark = kFALSE;

      TBranch *digitsBranch = 0;
      TTree *digitsTree;

      // Create an empty SPD digit array
      TObjArray *digitsArraySPD = new TObjArray();

      // Get the digits in update mode (we want to modify them if there are noisy channels)
      fITSLoader->UnloadDigits();
      fITSLoader->LoadDigits("update");

      // Get the number of events
      UInt_t numberOfEvents = (fRunLoader->TreeE()) ? static_cast<UInt_t>(fRunLoader->TreeE()->GetEntries()) : 0;

      // Loop over all events
      for (UInt_t event = 0; event < numberOfEvents; event++)
        {
          if (event%100 == 0) AliInfo(Form("Event #%d", event));

          // Get the current event
          fRunLoader->GetEvent(event);

          // Get the ITS digits tree
          digitsTree = fITSLoader->TreeD();

          // Get SPD branch that will contain all digits with marked noisy channels
          digitsBranch = digitsTree->GetBranch("ITSDigitsSPD");
          digitsBranch->SetAddress(&digitsArraySPD);

          // Get the stored number of modules
          UInt_t numberOfModules = (Int_t)digitsTree->GetEntries();
          TObjArray **newDigitsArraySPD = new TObjArray*[numberOfModules];

          Int_t *digitNumber = new Int_t[numberOfModules];
          for (UInt_t m = 0; m < numberOfModules; m++)
            {
              newDigitsArraySPD[m] = new TObjArray();
              digitNumber[m] = 0;
            }

          AliInfo(Form("ent = %d", (Int_t)digitsTree->GetEntries()));

          // Reset the SPD digit arrays to make sure they are empty
          digitsArraySPD->Clear();

          // Get the SPD digits branch from the original digits tree and set the address
          digitsBranch = digitsTree->GetBranch("ITSDigitsSPD");
          digitsBranch->SetAddress(&digitsArraySPD);

          // Loop over all modules
          for (UInt_t module = 0; module < fNumberOfModules; module++)
            {
              // Get event data for current module
              digitsTree->GetEvent(module);

              // Get the hits in the current module
              TObjArray *moduleObjArray = fBadChannelsContainer->CreateModuleObjArray(module);

              // Get the number of entries
              numberOfDigits = digitsArraySPD->GetEntries();

              // Loop over all digits and all channels
              for (UInt_t digit = 0; digit < numberOfDigits; digit++)
                {
                  // Get the current digit
                  digitSPD = (AliITSdigitSPD*) digitsArraySPD->At(digit);
                  newDigit[0] = digitSPD->GetCoord1(); // row
                  newDigit[1] = digitSPD->GetCoord2(); // column
                  newDigit[2] = digitSPD->GetSignal(); // signal

                  // Check if this channel is noisy
                  // (Compare with all stored channels in the bad channels array)
                  if (aux->Find(digitSPD, moduleObjArray))
                    {
                      // Set the mark flag and break the loop
                      mark = kTRUE;
                    }

                  // Store this digit in the SPD digits array using a placement new operation
                  new ((*newDigitsArraySPD[module])[digitNumber[module]]) AliITSdigitSPD(newDigit);

                  // Mark it if noisy and store in the noisy channel array
                  if (mark)
                    {
                      // Store this digit in the marked SPD digits array using a placement new operation
                      //new ((*badChannels[m])[numberOfBadChannels[m]]) AliITSChannelSPD(newBadChannel);
                      //new ((*newDigitsArraySPD[module])[digitNumber[module]]) AliITSdigitSPD(newDigit);

                      // Mark the original channel as noisy
                      ((*newDigitsArraySPD[module])[digitNumber[module]])->SetBit(kNoisyChannel);

                      mark = kFALSE;
                    }

                  digitNumber[module]++;

                } // end digit loop

              // Cleanup
              delete moduleObjArray;
              moduleObjArray = 0;

            } // end module loop

          digitsBranch->Reset();
          digitsBranch->ResetAddress();

          // Cleanup
          delete digitsArraySPD;
          digitsArraySPD = 0;
          digitsTree->Reset();

          // WHY THIS RANGE?????????????????????????????????????????????????????????????????????
          for (UInt_t n = 0; n < event; n++)
            {
              digitsTree->SetBranchAddress("ITSDigitsSPD", &newDigitsArraySPD[n]);
              digitsTree->Fill();
            }

          digitsTree->AutoSave();

          // Cleanup
          for (UInt_t n = 0; n < event; n++)
            {
              delete newDigitsArraySPD[n];
            }
          delete [] newDigitsArraySPD;
          newDigitsArraySPD = 0;
          delete [] digitNumber;
          digitNumber = 0;
          delete digitsTree;
          digitsTree = 0;

        } // end loop over all events

      // Unload the digits
      fITSLoader->UnloadDigits();

      // Cleanup
      delete aux;
      aux = 0;
    }
}

/*

//__________________________________________________________________________
Bool_t AliITSPreprocessorSPD::Store(AliCDBId &id, AliCDBMetaData *md)
{
  // Store the bad channels object in the calibration database
  // (See the corresponding run macro for further explanations)
  //
  // Input : fBadChannelsObjArray (now containing all found bad channels), object meta data
  // Output: Database file containing the bad channels
  // Return: kTRUE if successful

  Bool_t status = kFALSE;

  AliInfo("Storing bad channels");

  // Add the bad channels array to the storage container
  // (ownership is passed to the AliRunDataStorage object)
  fBadChannelsContainer = new AliITSBadChannelsSPD();

  // Convert the bad channels from TObjArray to Int_t[]
  AliITSPreprocessorSPD::ConvertObjToIntArray();

  // Store the arrays in the bad channels container object
  const Int_t kBadChannelsArraySize =
    2*fBadChannelsObjArray->GetEntries() + fNumberOfModules;
  fBadChannelsContainer->Put(fBadChannelsIntArray, kBadChannelsArraySize,
                             fBadChannelsIndexArray, fNumberOfModules);

  // Store the container
 if(!AliCDBManager::Instance()->IsDefaultStorageSet()) {
   //AliError("No storage set!");
   //   return status;
   AliCDBManager::Instance()->SetDefaultStorage("local://Calib");
 }
  
  if (AliCDBManager::Instance()->GetDefaultStorage()->Put(fBadChannelsContainer, id, md))
    {
      AliInfo("Bad channels object stored in database");
      status = kTRUE;
    }
  else
    {
      AliError("Failed to store object in database");
    }

  return status;
}

*/
//__________________________________________________________________________
void AliITSPreprocessorSPD::ConvertObjToIntArray()
{
  // Convert the bad channel TObjArray to an Int_t array
  //
  // Input : fBadChannelsObjArray (now containing all found bad channels)
  // Output: fBadChannelsIntArray
  // Return: (void)
  //
  // Data encoding:
  //                The TObjArray of this class (fBadChannelsObjArray) is converted to a sequential
  //                Int_t array (fBadChannelsIntArray) in this method. For each module, the first
  //                stored number is the number of bad channels in the current module. This is
  //                followed by all the columns and rows of the bad channels:
  //
  //                badChannelsArray =
  //                | N(m) | col0 | row0 | .. | colN(m) | N(m+1) | col0 | row0 | ...
  //                .   .......... module m .........   .   .... module m+1 ......
  //
  //                The bad channels index array (fBadChannelsIndexArray) contains the indices of
  //                the badChannelsArray, i.e. where the bad channels in certain module starts:
  //
  //                fBadChannelsObjArray =
  //                | i0 | i1 | .. | iM | (where M = the number of SPD modules)
  //
  //                e.g. i1 corresponds to the index of the badChannelsArray where N(1) is stored,
  //                i.e. the number of bad channels for module 1

  const Int_t kBadChannelsArraySize =
    2*fBadChannelsObjArray->GetEntries() + fNumberOfModules;
  fBadChannelsIntArray = new Int_t[kBadChannelsArraySize]; // Will be deleted in dtor
  AliITSChannelSPD *channel = 0;
  Int_t i = 0;
  Int_t j = 0;
  Int_t k = 0;

  // Loop over all modules
  for (UInt_t module = 0; module < fNumberOfModules; module++)
    {
      // Encode the number of bad channels of the current module
      fBadChannelsIntArray[k++] = fNumberOfBadChannels[module];

      // The columns and rows of the fBadChannelsObjArray will be stored sequentially
      // in the Int_t array
      j = 0;
      while (j < fNumberOfBadChannels[module])
        {
          channel = (AliITSChannelSPD *) fBadChannelsObjArray->At(i++);
          fBadChannelsIntArray[k++] = channel->GetColumn();
          fBadChannelsIntArray[k++] = channel->GetRow();

          // Go to next bad channel
          j++;
        }
    }
}


//__________________________________________________________________________
Bool_t AliITSPreprocessorSPD::Store(AliCDBId& /*id*/, AliCDBMetaData* /*md*/, Int_t runNumber)
{
  // Store the bad channels object in the calibration database
  // (See the corresponding run macro for further explanations)
  //
  // Input : fBadChannelsObjArray (now containing all found bad channels), object meta data
  // Output: Database file containing the bad channels
  // Return: kTRUE if successful

  Bool_t status = kFALSE;

  AliInfo("Storing bad channels");
 
  if(!AliCDBManager::Instance()->IsDefaultStorageSet()) {
    AliWarning("No storage set! Will use dummy one");
    AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
  }

  
  AliCDBEntry *entrySPD = AliCDBManager::Instance()->Get("ITS/Calib/CalibSPD", runNumber);
  if(!entrySPD){
    AliWarning("Calibration object retrieval failed! Dummy calibration will be used.");
    AliCDBStorage *origStorage = AliCDBManager::Instance()->GetDefaultStorage();
    AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT");
        
    entrySPD = AliCDBManager::Instance()->Get("ITS/Calib/CalibSPD", runNumber);
    AliCDBManager::Instance()->SetDefaultStorage(origStorage);
  }

  TObjArray *respSPD = (TObjArray *)entrySPD->GetObject();

  if ((! respSPD)) {
    AliWarning("Can not get calibration from calibration database !");
    return kFALSE;
  }
  
  Int_t i=0;
  AliITSChannelSPD *channel = 0;
  AliITSCalibrationSPD* res;
  for (Int_t module=0; module<respSPD->GetEntries(); module++) {
    Int_t k=0;
    res = (AliITSCalibrationSPD*) respSPD->At(module);
    res->SetNBadChannels(fNumberOfBadChannels[module]*2+1);
    res->AddBadChannel(fNumberOfBadChannels[module],k++);
    Int_t j = 0;
    while (j < fNumberOfBadChannels[module])
        {
          channel = (AliITSChannelSPD *) fBadChannelsObjArray->At(i++);
          res->AddBadChannel(channel->GetColumn(),k++);
          res->AddBadChannel(channel->GetRow(),k++);

          // Go to next bad channel
          j++;
        }
    
    
  }
  
  AliCDBManager::Instance()->Put(entrySPD);
  entrySPD->SetObject(NULL);
  entrySPD->SetOwner(kTRUE);

  delete entrySPD;
  status=kTRUE;
  return status;
}