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