EMCAL/AliCaloCalibSignal.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 /* $Id: AliCaloCalibSignal.cxx $ */
  16
  17 //________________________________________________________________________
  18 //
  19 // A help class for monitoring and calibration tools: MOOD, AMORE etc.,
  20 // It can be created and used a la (ctor):
  21 /*
  22   //Create the object for making the histograms
  23   fSignals = new AliCaloCalibSignal( fDetType );
  24   // AliCaloCalibSignal knows how many modules we have for PHOS or EMCAL
  25   fNumModules = fSignals->GetModules();
  26 */
  27 // fed an event:
  28 //  fSignals->ProcessEvent(fCaloRawStream,fRawEventHeaderBase);
  29 // asked to draw graphs or profiles:
  30 //  fSignals->GetGraphAmpVsTimeHighGain(module,column,row)->Draw("ap");
  31 // or
  32 //  fSignals->GetProfAmpVsTimeHighGain(module,column,row)->Draw();
  33 // etc.
  34 //________________________________________________________________________
  35
  36 #include "TFile.h"
  37
  38 #include "AliRawReader.h"
  39 #include "AliRawEventHeaderBase.h"
  40 #include "AliCaloRawStream.h"
  41
  42 //The include file
  43 #include "AliCaloCalibSignal.h"
  44
  45 ClassImp(AliCaloCalibSignal)
  46
  47 using namespace std;
  48
  49 // ctor; initialize everything in order to avoid compiler warnings
  50 // put some reasonable defaults
  51 AliCaloCalibSignal::AliCaloCalibSignal(kDetType detectorType) :
  52   TObject(),
  53   fDetType(kNone),
  54   fColumns(0),
  55   fRows(0),
  56   fModules(0),
  57   fCaloString(),
  58   fMapping(NULL),
  59   fRunNumber(-1),
  60   fStartTime(0),
  61   fAmpCut(50),
  62   fReqFractionAboveAmpCutVal(0.8),
  63   fReqFractionAboveAmp(kTRUE),
  64   fHour(0),
  65   fLatestHour(0),
  66   fUseAverage(kTRUE),
  67   fSecInAverage(1800),
  68   fNEvents(0),
  69   fNAcceptedEvents(0)
  70 {
  71   //Default constructor. First we set the detector-type related constants.
  72   if (detectorType == kPhos) {
  73     fColumns = fgkPhosCols;
  74     fRows = fgkPhosRows;
  75     fModules = fgkPhosModules;
  76     fCaloString = "PHOS";
  77   }
  78   else {
  79     //We'll just trust the enum to keep everything in line, so that if detectorType
  80     //isn't kPhos then it is kEmCal. Note, however, that this is not necessarily the
  81     //case, if someone intentionally gives another number
  82     fColumns = fgkEmCalCols;
  83     fRows = fgkEmCalRows;
  84     fModules = fgkEmCalModules;
  85     fCaloString = "EMCAL";
  86   }
  87
  88   fDetType = detectorType;
  89
  90   // Set the number of points for each Amp vs. Time graph to 0
  91   memset(fNHighGain, 0, sizeof(fNHighGain));
  92   memset(fNLowGain, 0, sizeof(fNLowGain));
  93
  94   CreateGraphs(); // set up the TGraphs
  95
  96   // init TProfiles to NULL=0 also
  97   memset(fProfAmpVsTimeHighGain, 0, sizeof(fProfAmpVsTimeHighGain));
  98   memset(fProfAmpVsTimeLowGain, 0, sizeof(fProfAmpVsTimeLowGain));
  99 }
 100
 101 // dtor
 102 //_____________________________________________________________________
 103 AliCaloCalibSignal::~AliCaloCalibSignal()
 104 {
 105   ClearObjects();
 106 }
 107
 108 //_____________________________________________________________________
 109 void AliCaloCalibSignal::ClearObjects()
 110 {
 111   // delete what was created in the ctor (TGraphs), and possible later (TProfiles)
 112   for (int i=0; i<fgkMaxTowers; i++) {
 113     if ( fGraphAmpVsTimeHighGain[i] ) { delete fGraphAmpVsTimeHighGain[i]; }
 114     if ( fGraphAmpVsTimeLowGain[i] ) { delete fGraphAmpVsTimeLowGain[i]; }
 115     if ( fProfAmpVsTimeHighGain[i] ) { delete fProfAmpVsTimeHighGain[i]; }
 116     if ( fProfAmpVsTimeLowGain[i] ) { delete fProfAmpVsTimeLowGain[i]; }
 117   }
 118   // set pointers
 119   memset(fGraphAmpVsTimeHighGain, 0, sizeof(fGraphAmpVsTimeHighGain));
 120   memset(fGraphAmpVsTimeLowGain, 0, sizeof(fGraphAmpVsTimeLowGain));
 121   memset(fProfAmpVsTimeHighGain, 0, sizeof(fProfAmpVsTimeHighGain));
 122   memset(fProfAmpVsTimeLowGain, 0, sizeof(fProfAmpVsTimeLowGain));
 123
 124   return;
 125 }
 126
 127 // copy ctor
 128 //_____________________________________________________________________
 129 AliCaloCalibSignal::AliCaloCalibSignal(const AliCaloCalibSignal &sig) :
 130   TObject(sig),
 131   fDetType(sig.GetDetectorType()),
 132   fColumns(sig.GetColumns()),
 133   fRows(sig.GetRows()),
 134   fModules(sig.GetModules()),
 135   fCaloString(sig.GetCaloString()),
 136   fMapping(NULL), //! note that we are not copying the map info
 137   fRunNumber(sig.GetRunNumber()),
 138   fStartTime(sig.GetStartTime()),
 139   fAmpCut(sig.GetAmpCut()),
 140   fReqFractionAboveAmpCutVal(sig.GetReqFractionAboveAmpCutVal()),
 141   fReqFractionAboveAmp(sig.GetReqFractionAboveAmp()),
 142   fHour(sig.GetHour()),
 143   fLatestHour(sig.GetLatestHour()),
 144   fUseAverage(sig.GetUseAverage()),
 145   fSecInAverage(sig.GetSecInAverage()),
 146   fNEvents(sig.GetNEvents()),
 147   fNAcceptedEvents(sig.GetNAcceptedEvents())
 148 {
 149   // also the TGraph contents
 150   AddInfo(&sig);
 151 }
 152
 153 // assignment operator; use copy ctor to make life easy..
 154 //_____________________________________________________________________
 155 AliCaloCalibSignal& AliCaloCalibSignal::operator = (const AliCaloCalibSignal &source)
 156 {
 157   // assignment operator; use copy ctor
 158   if (&source == this) return *this;
 159
 160   new (this) AliCaloCalibSignal(source);
 161   return *this;
 162 }
 163
 164 //_____________________________________________________________________
 165 void AliCaloCalibSignal::CreateGraphs()
 166 {
 167   //Then, loop for the requested number of modules
 168   TString title, name;
 169   for (int i = 0; i < fModules; i++) {
 170
 171     // Amplitude vs. Time graph for each channel
 172     for(int ic=0;ic < fColumns;ic++){
 173       for(int ir=0;ir < fRows;ir++){
 174
 175         int id = GetTowerNum(i, ic, ir);
 176
 177         // high gain graph
 178         name = "fGraphAmpVsTimeHighGain_"; name += i;
 179         name += "_"; name += ic;
 180         name += "_"; name += ir;
 181         title = "Amp vs. Time High Gain Mod "; title += i;
 182         title += " Col "; title += ic;
 183         title += " Row "; title += ir;
 184
 185         fGraphAmpVsTimeHighGain[id] = new TGraph();
 186         fGraphAmpVsTimeHighGain[id]->SetName(name);
 187         fGraphAmpVsTimeHighGain[id]->SetTitle(title);
 188         fGraphAmpVsTimeHighGain[id]->SetMarkerStyle(20);
 189
 190         // Low Gain
 191         name = "fGraphAmpVsTimeLowGain_"; name += i;
 192         name += "_"; name += ic;
 193         name += "_"; name += ir;
 194         title = "Amp vs. Time Low Gain Mod "; title += i;
 195         title += " Col "; title += ic;
 196         title += " Row "; title += ir;
 197
 198         fGraphAmpVsTimeLowGain[id] = new TGraph();
 199         fGraphAmpVsTimeLowGain[id]->SetName(name);
 200         fGraphAmpVsTimeLowGain[id]->SetTitle(title);
 201         fGraphAmpVsTimeLowGain[id]->SetMarkerStyle(20);
 202
 203       }
 204     }
 205
 206   }//end for nModules
 207 }
 208
 209 //_____________________________________________________________________
 210 void AliCaloCalibSignal::Reset()
 211 {
 212   Zero(); // set all counters to 0
 213   ClearObjects(); // delete previous TGraphs and TProfiles
 214   CreateGraphs(); // and create some new ones
 215   return;
 216 }
 217
 218 //_____________________________________________________________________
 219 void AliCaloCalibSignal::Zero()
 220 {
 221   // set all counters to 0; not cuts etc.though
 222   fHour = 0;
 223   fLatestHour = 0;
 224   fNEvents = 0;
 225   fNAcceptedEvents = 0;
 226   return;
 227 }
 228
 229 //_____________________________________________________________________
 230 Bool_t AliCaloCalibSignal::CheckFractionAboveAmp(int *AmpVal, int nTotChan)
 231 {
 232   int nAbove = 0;
 233
 234   int TowerNum = 0;
 235   for (int i = 0; i<fModules; i++) {
 236     for (int j = 0; j<fColumns; j++) {
 237       for (int k = 0; k<fRows; k++) {
 238         TowerNum = GetTowerNum(i,j,k);
 239         if (AmpVal[TowerNum] > fAmpCut) {
 240           nAbove++;
 241         }
 242       }
 243     }
 244   }
 245
 246   double fraction = (1.0*nAbove) / nTotChan;
 247
 248   if (fraction > fReqFractionAboveAmpCutVal) {
 249     return true;
 250   }
 251   else return false;
 252 }
 253
 254 //_____________________________________________________________________
 255 Bool_t AliCaloCalibSignal::AddInfo(const AliCaloCalibSignal *sig)
 256 {
 257   // just do this for the basic graphs/profiles that get filled in ProcessEvent
 258   // may not have data for all channels, but let's just Add everything..
 259   // Note: this method will run into problems with TProfile adding if the binning of
 260   // the local profiles is not the same as those provided by the argument *sig..
 261   int numGraphPoints = 0;
 262   int id = 0;
 263   int ip = 0;
 264   for (int i = 0; i < fModules; i++) {
 265     for (int j = 0; j < fColumns; j++) {
 266       for (int k = 0; k < fRows; k++) {
 267
 268         id = GetTowerNum(i,j,k);
 269
 270         if(fUseAverage){ // add to Profiles
 271           if (sig->GetProfAmpVsTimeHighGain(id)) {
 272             GetProfAmpVsTimeHighGain(id)->Add(sig->GetProfAmpVsTimeHighGain(id));
 273           }
 274           if (sig->GetProfAmpVsTimeLowGain(id)) {
 275             GetProfAmpVsTimeLowGain(id)->Add(sig->GetProfAmpVsTimeLowGain(id));
 276           }
 277         }
 278         else{ // add to Graphs
 279           // high gain
 280           numGraphPoints= sig->GetGraphAmpVsTimeHighGain(id)->GetN();
 281           if (numGraphPoints > 0) {
 282             // get the values
 283             double *graphX = sig->GetGraphAmpVsTimeHighGain(id)->GetX();
 284             double *graphY = sig->GetGraphAmpVsTimeHighGain(id)->GetY();
 285             for(ip=0; ip < numGraphPoints; ip++){
 286               fGraphAmpVsTimeHighGain[id]->SetPoint(fNHighGain[id]++,graphX[ip],graphY[ip]);
 287             }
 288           }
 289           // low gain
 290           numGraphPoints= sig->GetGraphAmpVsTimeLowGain(id)->GetN();
 291           if (numGraphPoints > 0) {
 292             // get the values
 293             double *graphX = sig->GetGraphAmpVsTimeLowGain(id)->GetX();
 294             double *graphY = sig->GetGraphAmpVsTimeLowGain(id)->GetY();
 295             for(ip=0; ip < numGraphPoints; ip++){
 296               fGraphAmpVsTimeLowGain[id]->SetPoint(fNLowGain[id]++,graphX[ip],graphY[ip]);
 297             }
 298           }
 299
 300         }
 301
 302       }//end for nModules
 303     }//end for nColumns
 304   }//end for nRows
 305
 306   return kTRUE;//We succesfully added info from the supplied object
 307 }
 308
 309 //_____________________________________________________________________
 310 Bool_t AliCaloCalibSignal::ProcessEvent(AliRawReader *rawReader)
 311 {
 312   // if fMapping is NULL the rawstream will crate its own mapping
 313   AliCaloRawStream rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
 314
 315   return ProcessEvent( &rawStream, (AliRawEventHeaderBase*)rawReader->GetEventHeader() );
 316 }
 317
 318 //_____________________________________________________________________
 319 Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStream *in, AliRawEventHeaderBase *aliHeader)
 320 {
 321   // Method to process=analyze one event in the data stream
 322   if (!in) return kFALSE; //Return right away if there's a null pointer
 323
 324   fNEvents++; // one more event
 325
 326   // PHOS has more towers than EMCAL, so use PHOS numbers to set array sizes
 327   int AmpValHighGain[fgkMaxTowers];
 328   int AmpValLowGain[fgkMaxTowers];
 329
 330   memset(AmpValHighGain, 0, sizeof(AmpValHighGain));
 331   memset(AmpValLowGain, 0, sizeof(AmpValLowGain));
 332
 333   int sample, i = 0; //The sample temp, and the sample number in current event.
 334   int max = fgkSampleMin, min = fgkSampleMax;//Use these for picking the signal
 335   int gain = 0;
 336
 337   // Number of Low and High gain channels for this event:
 338   int nLowChan = 0;
 339   int nHighChan = 0;
 340
 341   int TowerNum = 0; // array index for TGraphs etc.
 342
 343   // loop first to get the fraction of channels with amplitudes above cut
 344   while (in->Next()) {
 345     sample = in->GetSignal(); //Get the adc signal
 346     if (sample < min) min = sample;
 347     if (sample > max) max = sample;
 348     i++;
 349     if ( i >= in->GetTimeLength()) {
 350       //If we're here then we're done with this tower
 351       gain = 1 - in->IsLowGain();
 352
 353       int arrayPos = in->GetModule(); //The modules are numbered starting from 0
 354       if (arrayPos >= fModules) {
 355         //TODO: return an error message, if appopriate (perhaps if debug>0?)
 356         return kFALSE;
 357       }
 358
 359       //Debug
 360       if (arrayPos < 0 || arrayPos >= fModules) {
 361         printf("Oh no: arrayPos = %i.\n", arrayPos);
 362       }
 363
 364       // get tower number for AmpVal array
 365       TowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow());
 366
 367       if (gain == 0) {
 368         // fill amplitude into the array
 369         AmpValLowGain[TowerNum] = max - min;
 370         nLowChan++;
 371       }
 372       else if (gain==1) {//fill the high gain ones
 373         // fill amplitude into the array
 374         AmpValHighGain[TowerNum] = max - min;
 375         nHighChan++;
 376       }//end if gain
 377
 378
 379       max = fgkSampleMin; min = fgkSampleMax;
 380       i = 0;
 381
 382     }//End if end of tower
 383
 384   }//end while, of stream
 385
 386   // now check if it was a led event, only use high gain (that should be sufficient)
 387   if (fReqFractionAboveAmp) {
 388     bool ok = false;
 389     if (nHighChan > 0) {
 390       ok = CheckFractionAboveAmp(AmpValHighGain, nHighChan);
 391     }
 392     if (!ok) return false; // skip event
 393   }
 394
 395   fNAcceptedEvents++; // one more event accepted
 396
 397   if (fStartTime == 0) { // if start-timestamp wasn't set,we'll pick it up from the first event we encounter
 398     fStartTime = aliHeader->Get("Timestamp");
 399   }
 400
 401   fHour = (aliHeader->Get("Timestamp")-fStartTime)/(double)fgkNumSecInHr;
 402   if (fLatestHour < fHour) {
 403     fLatestHour = fHour;
 404   }
 405
 406   // it is a led event, now fill graphs (maybe profiles later)
 407   for(int i=0;i<fModules;i++){
 408     for(int j=0;j<fColumns;j++){
 409       for(int k=0;k<fRows;k++){
 410
 411         TowerNum = GetTowerNum(i, j, k);
 412
 413         if(AmpValHighGain[TowerNum]) {
 414           fGraphAmpVsTimeHighGain[TowerNum]->SetPoint(fNHighGain[TowerNum]++,fHour,AmpValHighGain[TowerNum]);
 415         }
 416         if(AmpValLowGain[TowerNum]) {
 417           fGraphAmpVsTimeLowGain[TowerNum]->SetPoint(fNLowGain[TowerNum]++,fHour,AmpValLowGain[TowerNum]);
 418         }
 419       }
 420     }
 421   }
 422
 423   return kTRUE;
 424 }
 425
 426 //_____________________________________________________________________
 427 void AliCaloCalibSignal::CreateProfile(int imod, int ic, int ir, int towerId, int gain,
 428                                        int nbins, double min, double max)
 429 { //! create/setup a TProfile
 430   TString title, name;
 431   if (gain == 0) {
 432     name = "fProfAmpVsTimeLowGain_";
 433     title = "Amp vs. Time Low Gain Mod ";
 434   }
 435   else if (gain == 1) {
 436     name = "fProfAmpVsTimeHighGain_";
 437     title = "Amp vs. Time High Gain Mod ";
 438   }
 439   name += imod;
 440   name += "_"; name += ic;
 441   name += "_"; name += ir;
 442   title += imod;
 443   title += " Col "; title += ic;
 444   title += " Row "; title += ir;
 445
 446   // use "s" option for RMS
 447   if (gain==0) {
 448     fProfAmpVsTimeLowGain[towerId] = new TProfile(name,title, nbins, min, max,"s");
 449   }
 450   else if (gain==1) {
 451     fProfAmpVsTimeHighGain[towerId] = new TProfile(name,title, nbins, min, max,"s");
 452   }
 453
 454   return;
 455 }
 456 //_____________________________________________________________________
 457 Bool_t AliCaloCalibSignal::Save(TString fileName, Bool_t saveEmptyGraphs)
 458 {
 459   //Saves all the histograms (or profiles, to be accurate) to the designated file
 460
 461   TFile destFile(fileName, "recreate");
 462
 463   if (destFile.IsZombie()) {
 464     return kFALSE;
 465   }
 466
 467   destFile.cd();
 468
 469   // setup variables for the TProfile plot
 470   int numProfBins = 0;
 471   double timeMin = 0;
 472   double timeMax = 0;
 473   if (fUseAverage) {
 474     if (fSecInAverage > 0) {
 475       numProfBins = (int)( (fLatestHour*fgkNumSecInHr)/fSecInAverage + 1 ); // round-off
 476     }
 477     numProfBins += 2; // add extra buffer : first and last
 478     double binSize = 1.0*fSecInAverage / fgkNumSecInHr;
 479     timeMin = - binSize;
 480     timeMax = timeMin + numProfBins*binSize;
 481   }
 482
 483   int numGraphPoints= 0;
 484   int TowerNum = 0;
 485   for (int i = 0; i < fModules; i++) {
 486
 487     for(int ic=0;ic < fColumns;ic++){
 488       for(int ir=0;ir < fRows;ir++){
 489
 490         TowerNum = GetTowerNum(i, ic, ir);
 491
 492         // 1st: high gain
 493         numGraphPoints= fGraphAmpVsTimeHighGain[TowerNum]->GetN();
 494         if( numGraphPoints>0 || saveEmptyGraphs) {
 495
 496           // average the graphs points over time if requested and put them in a profile plot
 497           if(fUseAverage && numGraphPoints>0) {
 498
 499             // get the values
 500             double *graphX = fGraphAmpVsTimeHighGain[TowerNum]->GetX();
 501             double *graphY = fGraphAmpVsTimeHighGain[TowerNum]->GetY();
 502
 503             // create the TProfile: 1 is for High gain
 504             CreateProfile(i, ic, ir, TowerNum, 1,
 505                           numProfBins, timeMin, timeMax);
 506
 507             // loop over graph points and fill profile
 508             for(int ip=0; ip < numGraphPoints; ip++){
 509               fProfAmpVsTimeHighGain[TowerNum]->Fill(graphX[ip],graphY[ip]);
 510             }
 511
 512             fProfAmpVsTimeHighGain[TowerNum]->GetXaxis()->SetTitle("Hours");
 513             fProfAmpVsTimeHighGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal");
 514             fProfAmpVsTimeHighGain[TowerNum]->Write();
 515
 516           }
 517            else{
 518              //otherwise, just save the graphs and forget the profiling
 519              fGraphAmpVsTimeHighGain[TowerNum]->GetXaxis()->SetTitle("Hours");
 520              fGraphAmpVsTimeHighGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal");
 521              fGraphAmpVsTimeHighGain[TowerNum]->Write();
 522            }
 523
 524         } // low gain graph info should be saved in one form or another
 525
 526         // 2nd: now go to the low gain case
 527         numGraphPoints= fGraphAmpVsTimeLowGain[TowerNum]->GetN();
 528         if( numGraphPoints>0 || saveEmptyGraphs) {
 529
 530           // average the graphs points over time if requested and put them in a profile plot
 531           if(fUseAverage && numGraphPoints>0) {
 532
 533             double *graphX = fGraphAmpVsTimeLowGain[TowerNum]->GetX();
 534             double *graphY = fGraphAmpVsTimeLowGain[TowerNum]->GetY();
 535
 536             // create the TProfile: 0 is for Low gain
 537             CreateProfile(i, ic, ir, TowerNum, 0,
 538                           numProfBins, timeMin, timeMax);
 539
 540             // loop over graph points and fill profile
 541             for(int ip=0; ip < numGraphPoints; ip++){
 542               fProfAmpVsTimeLowGain[TowerNum]->Fill(graphX[ip],graphY[ip]);
 543             }
 544
 545             fProfAmpVsTimeLowGain[TowerNum]->GetXaxis()->SetTitle("Hours");
 546             fProfAmpVsTimeLowGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal");
 547             fProfAmpVsTimeLowGain[TowerNum]->Write();
 548
 549           }
 550           else{
 551              //otherwise, just save the graphs and forget the profiling
 552             fGraphAmpVsTimeLowGain[TowerNum]->GetXaxis()->SetTitle("Hours");
 553             fGraphAmpVsTimeLowGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal");
 554             fGraphAmpVsTimeLowGain[TowerNum]->Write();
 555           }
 556
 557         } // low gain graph info should be saved in one form or another
 558
 559       } // fRows
 560     } // fColumns
 561
 562   } // fModules
 563   destFile.Close();
 564
 565   return kTRUE;
 566 }