ITS/AliITSPreprocessorSPD.cxx

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