EMCAL/AliCaloCalibPedestal.cxx

   1 #include "TH1.h"
   2 #include "TFile.h"
   3
   4 #include <fstream>
   5 #include <iostream>
   6
   7 //The include file
   8 #include "AliCaloCalibPedestal.h"
   9
  10 ClassImp(AliCaloCalibPedestal)
  11
  12 using namespace std;
  13
  14 // ctor; initialize everything in order to avoid compiler warnings
  15 AliCaloCalibPedestal::AliCaloCalibPedestal(kDetType detectorType) :
  16   TObject(),
  17   fPedestalLowGain(),
  18   fPedestalHighGain(),
  19   fPeakMinusPedLowGain(),
  20   fPeakMinusPedHighGain(),
  21   fPedestalLowGainDiff(),
  22   fPedestalHighGainDiff(),
  23   fPeakMinusPedLowGainDiff(),
  24   fPeakMinusPedHighGainDiff(),
  25   fPedestalLowGainRatio(),
  26   fPedestalHighGainRatio(),
  27   fPeakMinusPedLowGainRatio(),
  28   fPeakMinusPedHighGainRatio(),
  29   fDeadMap(),
  30   fDeadTowers(0),
  31   fNewDeadTowers(0),
  32   fResurrectedTowers(0),
  33   fReference(0),
  34   fDetType(kNone),
  35   fColumns(0),
  36   fRows(0),
  37   fModules(0),
  38   fRunNumber(-1)
  39 {
  40   //Default constructor. First we set the detector-type related constants.
  41   if (detectorType == kPhos) {
  42     fColumns = fgkPhosCols;
  43     fRows = fgkPhosRows;
  44     fModules = fgkPhosModules;
  45   }
  46   else {
  47     //We'll just trust the enum to keep everything in line, so that if detectorType
  48     //isn't kPhos then it is kEmCal. Note, however, that this is not necessarily the
  49     //case, if someone intentionally gives another number
  50     fColumns = fgkEmCalCols;
  51     fRows = fgkEmCalRows;
  52     fModules = fgkEmCalModules;
  53   }
  54   fDetType = detectorType;
  55
  56   //Then, loop for the requested number of modules
  57   TString title, name;
  58   for (int i = 0; i < fModules; i++) {
  59     //Pedestals, low gain
  60     name = "hPedlowgain";
  61     name += i;
  62     title = "Pedestals, low gain, module ";
  63     title += i;
  64     fPedestalLowGain.Add(new TProfile2D(name, title,
  65                                         fColumns, 0.0, fColumns,
  66                                         fRows, -fRows, 0.0));
  67
  68     //Pedestals, high gain
  69     name = "hPedhighgain";
  70     name += i;
  71     title = "Pedestals, high gain, module ";
  72     title += i;
  73     fPedestalHighGain.Add(new TProfile2D(name, title,
  74                                          fColumns, 0.0, fColumns,
  75                                          fRows, -fRows, 0.0));
  76
  77     //Peak-Pedestals, low gain
  78     name = "hPeakMinusPedlowgain";
  79     name += i;
  80     title = "Peak-Pedestal, low gain, module ";
  81     title += i;
  82     fPeakMinusPedLowGain.Add(new TProfile2D(name, title,
  83                                             fColumns, 0.0, fColumns,
  84                                             fRows, -fRows, 0.0));
  85
  86     //Peak-Pedestals, high gain
  87     name = "hPeakMinusPedhighgain";
  88     name += i;
  89     title = "Peak-Pedestal, high gain, module ";
  90     title += i;
  91     fPeakMinusPedHighGain.Add(new TProfile2D(name, title,
  92                                              fColumns, 0.0, fColumns,
  93                                              fRows, -fRows, 0.0));
  94
  95     name = "hDeadMap";
  96     name += i;
  97     title = "Dead map, module ";
  98     title += i;
  99     fDeadMap.Add(new TH2D(name, title, fColumns, 0.0, fColumns,
 100                           fRows, -fRows, 0.0));
 101
 102   }//end for nModules create the histograms
 103
 104   //Compress the arrays, in order to remove the empty objects (a 16 slot array is created by default)
 105   fPedestalLowGain.Compress();
 106   fPedestalHighGain.Compress();
 107   fPeakMinusPedLowGain.Compress();
 108   fPeakMinusPedHighGain.Compress();
 109   fDeadMap.Compress();
 110   //Make them the owners of the profiles, so we don't need to care about deleting them
 111   //fPedestalLowGain.SetOwner();
 112   //fPedestalHighGain.SetOwner();
 113   //fPeakMinusPedLowGain.SetOwner();
 114   //fPeakMinusPedHighGain.SetOwner();
 115
 116 }
 117
 118 // dtor
 119 //_____________________________________________________________________
 120 AliCaloCalibPedestal::~AliCaloCalibPedestal()
 121 {
 122   if (fReference) delete fReference;//Delete the reference object, if it has been loaded
 123   //TObjArray will delete the histos/profiles when it is deleted.
 124 }
 125
 126 // copy ctor
 127 //_____________________________________________________________________
 128 AliCaloCalibPedestal::AliCaloCalibPedestal(const AliCaloCalibPedestal &ped) :
 129   TObject(ped),
 130   fPedestalLowGain(),
 131   fPedestalHighGain(),
 132   fPeakMinusPedLowGain(),
 133   fPeakMinusPedHighGain(),
 134   fPedestalLowGainDiff(),
 135   fPedestalHighGainDiff(),
 136   fPeakMinusPedLowGainDiff(),
 137   fPeakMinusPedHighGainDiff(),
 138   fPedestalLowGainRatio(),
 139   fPedestalHighGainRatio(),
 140   fPeakMinusPedLowGainRatio(),
 141   fPeakMinusPedHighGainRatio(),
 142   fDeadMap(),
 143   fDeadTowers(ped.GetDeadTowerCount()),
 144   fNewDeadTowers(ped.GetDeadTowerNew()),
 145   fResurrectedTowers(ped.GetDeadTowerResurrected()),
 146   fReference( 0 ), //! note that we do not try to copy the reference info here
 147   fDetType(ped.GetDetectorType()),
 148   fColumns(ped.GetColumns()),
 149   fRows(ped.GetRows()),
 150   fModules(ped.GetModules()),
 151   fRunNumber(ped.GetRunNumber())
 152 {
 153   // Then the ObjArray ones; we add the histograms rather than trying TObjArray = assignment
 154   //DS: this has not really been tested yet..
 155   for (int i = 0; i < fModules; i++) {
 156     fPedestalLowGain.Add( ped.GetPedProfileLowGain(i) );
 157     fPedestalHighGain.Add( ped.GetPedProfileHighGain(i) );
 158     fPeakMinusPedLowGain.Add( ped.GetPeakProfileLowGain(i) );
 159     fPeakMinusPedHighGain.Add( ped.GetPeakProfileHighGain(i) );
 160
 161     fDeadMap.Add( ped.GetDeadMap(i) );
 162   }//end for nModules
 163
 164   //Compress the arrays, in order to remove the empty objects (a 16 slot array is created by default)
 165   fPedestalLowGain.Compress();
 166   fPedestalHighGain.Compress();
 167   fPeakMinusPedLowGain.Compress();
 168   fPeakMinusPedHighGain.Compress();
 169   fDeadMap.Compress();
 170 }
 171
 172 // assignment operator; use copy ctor to make life easy..
 173 //_____________________________________________________________________
 174 AliCaloCalibPedestal& AliCaloCalibPedestal::operator = (const AliCaloCalibPedestal &source)
 175 {
 176   if (&source == this) return *this;
 177
 178   new (this) AliCaloCalibPedestal(source);
 179   return *this;
 180 }
 181
 182 //_____________________________________________________________________
 183 void AliCaloCalibPedestal::Reset()
 184 {
 185   for (int i = 0; i < fModules; i++) {
 186     GetPedProfileLowGain(i)->Reset();
 187     GetPedProfileHighGain(i)->Reset();
 188     GetPeakProfileLowGain(i)->Reset();
 189     GetPeakProfileHighGain(i)->Reset();
 190     GetDeadMap(i)->Reset();
 191
 192     if (!fPedestalLowGainDiff.IsEmpty()) {
 193       //This means that the comparison profiles have been created.
 194
 195       GetPedProfileLowGainDiff(i)->Reset();
 196       GetPedProfileHighGainDiff(i)->Reset();
 197       GetPeakProfileLowGainDiff(i)->Reset();
 198       GetPeakProfileHighGainDiff(i)->Reset();
 199
 200       GetPedProfileLowGainRatio(i)->Reset();
 201       GetPedProfileHighGainRatio(i)->Reset();
 202       GetPeakProfileLowGainRatio(i)->Reset();
 203       GetPeakProfileHighGainRatio(i)->Reset();
 204     }
 205   }
 206   fDeadTowers = 0;
 207   fNewDeadTowers = 0;
 208   fResurrectedTowers = 0;
 209
 210   //To think about: should fReference be deleted too?... let's not do it this time, at least...
 211 }
 212
 213 //_____________________________________________________________________
 214 Bool_t AliCaloCalibPedestal::ProcessEvent(AliCaloRawStream *in)
 215 {
 216   if (!in) return kFALSE; //Return right away if there's a null pointer
 217
 218   in->SetOldRCUFormat(kTRUE);
 219
 220   int sample, i = 0; //The sample temp, and the sample number in current event.
 221   int max = fgkSampleMin, min = fgkSampleMax;//Use these for picking the pedestal
 222   int gain = 0;
 223
 224   while (in->Next()) {
 225     sample = in->GetSignal(); //Get the adc signal
 226     if (sample < min) min = sample;
 227     if (sample > max) max = sample;
 228     i++;
 229     if ( i >= in->GetTimeLength()) {
 230       //If we're here then we're done with this tower
 231       gain = 1 - in->IsLowGain();
 232
 233       int arrayPos = in->GetModule(); //The modules are numbered starting from 0
 234       if (arrayPos >= fModules) {
 235         //TODO: return an error message, if appopriate (perhaps if debug>0?)
 236         return kFALSE;
 237       }
 238
 239       //Debug
 240       if (arrayPos < 0 || arrayPos >= fModules) {
 241         printf("Oh no: arrayPos = %i.\n", arrayPos);
 242       }
 243
 244       //NOTE: coordinates are (column, row) for the profiles
 245       if (gain == 0) {
 246         //fill the low gain histograms
 247         ((TProfile2D*)fPedestalLowGain[arrayPos])->Fill(in->GetColumn(), -in->GetRow() - 1, min);
 248         ((TProfile2D*)fPeakMinusPedLowGain[arrayPos])->Fill(in->GetColumn(), -in->GetRow() - 1, max - min);
 249       }
 250       else {//fill the high gain ones
 251         ((TProfile2D*)fPedestalHighGain[arrayPos])->Fill(in->GetColumn(), -in->GetRow() - 1, min);
 252         ((TProfile2D*)fPeakMinusPedHighGain[arrayPos])->Fill(in->GetColumn(), -in->GetRow() - 1, max - min);
 253       }//end if gain
 254
 255       max = fgkSampleMin; min = fgkSampleMax;
 256       i = 0;
 257
 258     }//End if end of tower
 259
 260   }//end while, of stream
 261
 262   return kTRUE;
 263 }
 264
 265 //_____________________________________________________________________
 266 Bool_t AliCaloCalibPedestal::SaveHistograms(TString fileName, Bool_t saveEmptyHistos)
 267 {
 268   //Saves all the histograms (or profiles, to be accurate) to the designated file
 269
 270   TFile destFile(fileName, "recreate");
 271
 272   if (destFile.IsZombie()) {
 273     return kFALSE;
 274   }
 275
 276   destFile.cd();
 277
 278   for (int i = 0; i < fModules; i++) {
 279     if( ((TProfile2D *)fPeakMinusPedLowGain[i])->GetEntries() || saveEmptyHistos) {
 280       fPeakMinusPedLowGain[i]->Write();
 281     }
 282     if( ((TProfile2D *)fPeakMinusPedHighGain[i])->GetEntries() || saveEmptyHistos) {
 283       fPeakMinusPedHighGain[i]->Write();
 284     }
 285     if( ((TProfile2D *)fPedestalLowGain[i])->GetEntries() || saveEmptyHistos) {
 286       fPedestalLowGain[i]->Write();
 287     }
 288     if( ((TProfile2D *)fPedestalHighGain[i])->GetEntries() || saveEmptyHistos) {
 289       fPedestalHighGain[i]->Write();
 290     }
 291   }
 292
 293   destFile.Close();
 294
 295   return kTRUE;
 296 }
 297
 298 //_____________________________________________________________________
 299 Bool_t AliCaloCalibPedestal::LoadReferenceCalib(TString fileName, TString objectName)
 300 {
 301
 302   //Make sure that the histograms created when loading the object are not destroyed as the file object is destroyed
 303   TH1::AddDirectory(kFALSE);
 304
 305   TFile *sourceFile = new TFile(fileName);
 306   if (sourceFile->IsZombie()) {
 307     return kFALSE;//We couldn't load the reference
 308   }
 309
 310   if (fReference) delete fReference;//Delete the reference object, if it already exists
 311   fReference = 0;
 312
 313   fReference = (AliCaloCalibPedestal*)sourceFile->Get(objectName);
 314
 315   if (!fReference || !(fReference->InheritsFrom(AliCaloCalibPedestal::Class())) || (fReference->GetDetectorType() != fDetType)) {
 316     if (fReference) delete fReference;//Delete the object, in case we had an object of the wrong type
 317     fReference = 0;
 318     return kFALSE;
 319   }
 320
 321   delete sourceFile;
 322
 323   //Reset the histogram ownership behaviour. NOTE: a better workaround would be good, since this may accidentally set AddDirectory to true, even
 324   //if we are called by someone who has set it to false...
 325   TH1::AddDirectory(kTRUE);
 326
 327   return kTRUE;//We succesfully loaded the object
 328 }
 329
 330 //_____________________________________________________________________
 331 void AliCaloCalibPedestal::ValidateComparisonProfiles()
 332 {
 333   if (!fPedestalLowGainDiff.IsEmpty()) return; //The profiles already exist. We just check one, because they're all created at
 334   //the same time
 335
 336
 337   //Then, loop for the requested number of modules
 338   TString title, name;
 339   for (int i = 0; i < fModules; i++) {
 340     //Pedestals, low gain
 341     name = "hPedlowgainDiff";
 342     name += i;
 343     title = "Pedestals difference, low gain, module ";
 344     title += i;
 345     fPedestalLowGainDiff.Add(new TProfile2D(name, title,
 346                                             fColumns, 0.0, fColumns,
 347                                             fRows, -fRows, 0.0));
 348
 349     //Pedestals, high gain
 350     name = "hPedhighgainDiff";
 351     name += i;
 352     title = "Pedestals difference, high gain, module ";
 353     title += i;
 354     fPedestalHighGainDiff.Add(new TProfile2D(name, title,
 355                                              fColumns, 0.0, fColumns,
 356                                              fRows, -fRows, 0.0));
 357
 358     //Peak-Pedestals, low gain
 359     name = "hPeakMinusPedlowgainDiff";
 360     name += i;
 361     title = "Peak-Pedestal difference, low gain, module ";
 362     title += i;
 363     fPeakMinusPedLowGainDiff.Add(new TProfile2D(name, title,
 364                                                 fColumns, 0.0, fColumns,
 365                                                 fRows, -fRows, 0.0));
 366
 367     //Peak-Pedestals, high gain
 368     name = "hPeakMinusPedhighgainDiff";
 369     name += i;
 370     title = "Peak-Pedestal difference, high gain, module ";
 371     title += i;
 372     fPeakMinusPedHighGainDiff.Add(new TProfile2D(name, title,
 373                                                  fColumns, 0.0, fColumns,
 374                                                  fRows, -fRows, 0.0));
 375
 376     //Pedestals, low gain
 377     name = "hPedlowgainRatio";
 378     name += i;
 379     title = "Pedestals ratio, low gain, module ";
 380     title += i;
 381     fPedestalLowGainRatio.Add(new TProfile2D(name, title,
 382                                              fColumns, 0.0, fColumns,
 383                                              fRows, -fRows, 0.0));
 384
 385     //Pedestals, high gain
 386     name = "hPedhighgainRatio";
 387     name += i;
 388     title = "Pedestals ratio, high gain, module ";
 389     title += i;
 390     fPedestalHighGainRatio.Add(new TProfile2D(name, title,
 391                                               fColumns, 0.0, fColumns,
 392                                               fRows, -fRows, 0.0));
 393
 394     //Peak-Pedestals, low gain
 395     name = "hPeakMinusPedlowgainRatio";
 396     name += i;
 397     title = "Peak-Pedestal ratio, low gain, module ";
 398     title += i;
 399     fPeakMinusPedLowGainRatio.Add(new TProfile2D(name, title,
 400                                                  fColumns, 0.0, fColumns,
 401                                                  fRows, -fRows, 0.0));
 402
 403     //Peak-Pedestals, high gain
 404     name = "hPeakMinusPedhighgainRatio";
 405     name += i;
 406     title = "Peak-Pedestal ratio, high gain, module ";
 407     title += i;
 408     fPeakMinusPedHighGainRatio.Add(new TProfile2D(name, title,
 409                                                   fColumns, 0.0, fColumns,
 410                                                   fRows, -fRows, 0.0));
 411
 412   }//end for nModules create the histograms
 413 }
 414
 415 //_____________________________________________________________________
 416 void AliCaloCalibPedestal::ComputeDiffAndRatio()
 417 {
 418
 419   ValidateComparisonProfiles();//Make sure the comparison histos exist
 420
 421   if (!fReference) {
 422     return;//Return if the reference object isn't loaded
 423   }
 424
 425   for (int i = 0; i < fModules; i++) {
 426     //Compute the ratio of the histograms
 427
 428     ((TProfile2D*)fPedestalLowGainRatio[i])->Divide(GetPedProfileLowGain(i), fReference->GetPedProfileLowGain(i), 1.0, 1.0);
 429     ((TProfile2D*)fPedestalHighGainRatio[i])->Divide(GetPedProfileHighGain(i), fReference->GetPedProfileHighGain(i), 1.0, 1.0);
 430     ((TProfile2D*)fPeakMinusPedLowGainRatio[i])->Divide(GetPeakProfileLowGain(i), fReference->GetPeakProfileLowGain(i), 1.0, 1.0);
 431     ((TProfile2D*)fPeakMinusPedHighGainRatio[i])->Divide(GetPeakProfileHighGain(i), fReference->GetPeakProfileHighGain(i), 1.0, 1.0);
 432
 433     //For computing the difference, we cannot simply do TProfile2D->Add(), because that subtracts the sum of all entries,
 434     //which means that the mean of the new profile will not be the difference of the means. So do it by hand:
 435     for (int j = 0; j <= fColumns; j++) {
 436       for (int k = 0; k <= fRows; k++) {
 437         int bin = ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->GetBin(j+1, k+1);//Note that we assume here that all histos have the same structure...
 438         double diff = fReference->GetPeakProfileHighGain(i)->GetBinContent(bin) - GetPeakProfileHighGain(i)->GetBinContent(bin);
 439         ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinContent(j+1, k+1, diff);
 440         ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinEntries(bin, 1);
 441
 442         diff = fReference->GetPeakProfileLowGain(i)->GetBinContent(bin) - GetPeakProfileLowGain(i)->GetBinContent(bin);
 443         ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinContent(j+1, k+1, diff);
 444         ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinEntries(bin, 1);
 445
 446         diff = fReference->GetPedProfileHighGain(i)->GetBinContent(bin) - GetPedProfileHighGain(i)->GetBinContent(bin);
 447         ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinContent(j+1, k+1, diff);
 448         ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinEntries(bin, 1);
 449
 450         diff = fReference->GetPedProfileLowGain(i)->GetBinContent(bin) - GetPedProfileLowGain(i)->GetBinContent(bin);
 451         ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinContent(j+1, k+1, diff);
 452         ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinEntries(bin, 1);
 453
 454       } // rows
 455     } // columns
 456
 457   } // modules
 458
 459 }
 460
 461 //_____________________________________________________________________
 462 void AliCaloCalibPedestal::ComputeDeadTowers(int threshold, const char * deadMapFile)
 463 {//Computes the number of dead towers etc etc into memory, after this you can call the GetDead... -functions
 464   int countTot = 0;
 465   int countNew = 0;
 466   int countRes = 0;
 467   ofstream * fout = 0;
 468   ofstream * diff = 0;
 469   char name[512];//Quite a long temp buffer, just in case the filename includes a path
 470
 471   if (deadMapFile) {
 472     snprintf(name, 512, "%s.txt", deadMapFile);
 473     fout = new ofstream(name);
 474     snprintf(name, 512, "%sdiff.txt", deadMapFile);
 475     diff = new ofstream(name);
 476     if (!fout->is_open()) {
 477       delete fout;
 478       fout = 0;//Set the pointer to empty if the file was not opened
 479     }
 480     if (!diff->is_open()) {
 481       delete diff;
 482       fout = 0;//Set the pointer to empty if the file was not opened
 483     }
 484   }
 485
 486   for (int i = 0; i < fModules; i++) {
 487     if (GetPeakProfileHighGain(i)->GetEntries() > 0) { //don't care about empty histos
 488       for (int j = 1; j <= fColumns; j++) {
 489         for (int k = 1; k <= fRows; k++) {
 490
 491           if (GetPeakProfileHighGain(i)->GetBinContent(j, k) < threshold) {//It's dead
 492             countTot++;//One more dead total
 493             if (fout) {
 494               (*fout) << i << " "
 495                       << (fRows - k) << " "
 496                       << (j-1) << " "
 497                       << "1" << " "
 498                       << "0" << endl;//Write the status to the deadmap file, if the file is open.
 499             }
 500
 501             if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) >= threshold) {
 502               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kRecentlyDeceased);
 503               countNew++;//This tower wasn't dead before!
 504               if (diff) {
 505                 ( *diff) << i << " "
 506                          << (fRows - k) << " "
 507                          << (j - 1) << " "
 508                          << "1" << " "
 509                          << "0" << endl;//Write the status to the deadmap difference file, if the file is open.
 510               }
 511             }
 512             else {
 513               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kDead);//This has been dead before. Nothing new
 514             }
 515           }
 516           else { //It's ALIVE!!
 517             //Don't bother with writing the live ones.
 518             //if (fout)
 519             //  (*fout) << i << " "
 520             //     << (fRows - k) << " "
 521             //     << (j - 1) << " "
 522             //     << "1" << " "
 523             //     << "1" << endl;//Write the status to the deadmap file, if the file is open.
 524             if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) < threshold) {
 525               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kResurrected);
 526               countRes++; //This tower was dead before => it's a miracle! :P
 527               if (diff) {
 528                 (*diff) << i << " "
 529                         << (fRows - k) << " "
 530                         << (j - 1) << " "
 531                         << "1" << " "
 532                         << "1" << endl;//Write the status to the deadmap difference file, if the file is open.
 533               }
 534             }
 535             else {
 536               ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kAlive);
 537             }
 538           }
 539
 540         }//end for k/rows
 541       }//end for j/columns
 542     }//end if GetEntries >= 0
 543
 544   }//end for modules
 545
 546  if (fout) {
 547    fout->close();
 548    delete fout;
 549  }
 550
 551  fDeadTowers = countTot;
 552  fNewDeadTowers = countNew;
 553  fResurrectedTowers = countRes;
 554 }
 555