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 // get some info:
  30 //  fSignals->GetXXX..()
  31 // etc.
  32 //________________________________________________________________________
  33 #include <string>
  34 #include <sstream>
  35 #include <fstream>
  36
  37 #include "TProfile.h"
  38 #include "TFile.h"
  39
  40 #include "AliRawReader.h"
  41 #include "AliCaloRawStreamV3.h"
  42
  43 //The include file
  44 #include "AliCaloCalibSignal.h"
  45
  46 ClassImp(AliCaloCalibSignal)
  47
  48 using namespace std;
  49
  50 // variables for TTree filling; not sure if they should be static or not
  51 static int fChannelNum = 0; // for regular towers
  52 static int fRefNum = 0; // for LED
  53 static double fAmp = 0;
  54 static double fAvgAmp = 0;
  55 static double fRMS = 0;
  56
  57 // ctor; initialize everything in order to avoid compiler warnings
  58 // put some reasonable defaults
  59 AliCaloCalibSignal::AliCaloCalibSignal(kDetType detectorType) :
  60   TObject(),
  61   fDetType(kNone),
  62   fColumns(0),
  63   fRows(0),
  64   fLEDRefs(0),
  65   fModules(0),
  66   fCaloString(),
  67   fMapping(NULL),
  68   fRunNumber(-1),
  69   fStartTime(0),
  70   fAmpCut(40), // min. 40 ADC counts as default
  71   fReqFractionAboveAmpCutVal(0.6), // 60% in a strip, per default
  72   fReqFractionAboveAmp(kTRUE),
  73   fAmpCutLEDRef(100), // min. 100 ADC counts as default
  74   fReqLEDRefAboveAmpCutVal(kTRUE),
  75   fHour(0),
  76   fLatestHour(0),
  77   fUseAverage(kTRUE),
  78   fSecInAverage(1800),
  79   fNEvents(0),
  80   fNAcceptedEvents(0),
  81   fTreeAmpVsTime(NULL),
  82   fTreeAvgAmpVsTime(NULL),
  83   fTreeLEDAmpVsTime(NULL),
  84   fTreeLEDAvgAmpVsTime(NULL)
  85 {
  86   //Default constructor. First we set the detector-type related constants.
  87   if (detectorType == kPhos) {
  88     fColumns = fgkPhosCols;
  89     fRows = fgkPhosRows;
  90     fLEDRefs = fgkPhosLEDRefs;
  91     fModules = fgkPhosModules;
  92     fCaloString = "PHOS";
  93   }
  94   else {
  95     //We'll just trust the enum to keep everything in line, so that if detectorType
  96     //isn't kPhos then it is kEmCal. Note, however, that this is not necessarily the
  97     //case, if someone intentionally gives another number
  98     fColumns = AliEMCALGeoParams::fgkEMCALCols;
  99     fRows = AliEMCALGeoParams::fgkEMCALRows;
 100     fLEDRefs = AliEMCALGeoParams::fgkEMCALLEDRefs;
 101     fModules = AliEMCALGeoParams::fgkEMCALModules;
 102     fCaloString = "EMCAL";
 103   }
 104
 105   fDetType = detectorType;
 106
 107   ResetInfo(); // trees and counters
 108 }
 109
 110 // dtor
 111 //_____________________________________________________________________
 112 AliCaloCalibSignal::~AliCaloCalibSignal()
 113 {
 114   DeleteTrees();
 115 }
 116
 117 //_____________________________________________________________________
 118 void AliCaloCalibSignal::DeleteTrees()
 119 {
 120   // delete what was created in the ctor (TTrees)
 121   if (fTreeAmpVsTime) delete fTreeAmpVsTime;
 122   if (fTreeAvgAmpVsTime) delete fTreeAvgAmpVsTime;
 123   if (fTreeLEDAmpVsTime) delete fTreeLEDAmpVsTime;
 124   if (fTreeLEDAvgAmpVsTime) delete fTreeLEDAvgAmpVsTime;
 125   // and reset pointers
 126   fTreeAmpVsTime = NULL;
 127   fTreeAvgAmpVsTime = NULL;
 128   fTreeLEDAmpVsTime = NULL;
 129   fTreeLEDAvgAmpVsTime = NULL;
 130
 131   return;
 132 }
 133
 134 // copy ctor
 135 //_____________________________________________________________________
 136 AliCaloCalibSignal::AliCaloCalibSignal(const AliCaloCalibSignal &sig) :
 137   TObject(sig),
 138   fDetType(sig.GetDetectorType()),
 139   fColumns(sig.GetColumns()),
 140   fRows(sig.GetRows()),
 141   fLEDRefs(sig.GetLEDRefs()),
 142   fModules(sig.GetModules()),
 143   fCaloString(sig.GetCaloString()),
 144   fMapping(), //! note that we are not copying the map info
 145   fRunNumber(sig.GetRunNumber()),
 146   fStartTime(sig.GetStartTime()),
 147   fAmpCut(sig.GetAmpCut()),
 148   fReqFractionAboveAmpCutVal(sig.GetReqFractionAboveAmpCutVal()),
 149   fReqFractionAboveAmp(sig.GetReqFractionAboveAmp()),
 150   fAmpCutLEDRef(sig.GetAmpCutLEDRef()),
 151   fReqLEDRefAboveAmpCutVal(sig.GetReqLEDRefAboveAmpCutVal()),
 152   fHour(sig.GetHour()),
 153   fLatestHour(sig.GetLatestHour()),
 154   fUseAverage(sig.GetUseAverage()),
 155   fSecInAverage(sig.GetSecInAverage()),
 156   fNEvents(sig.GetNEvents()),
 157   fNAcceptedEvents(sig.GetNAcceptedEvents()),
 158   fTreeAmpVsTime(),
 159   fTreeAvgAmpVsTime(),
 160   fTreeLEDAmpVsTime(),
 161   fTreeLEDAvgAmpVsTime()
 162 {
 163   // also the TTree contents
 164   AddInfo(&sig);
 165   for (Int_t i = 0; i<2*fgkMaxTowers; i++) {
 166     if(i < fgkMaxTowers){
 167       fNHighGain[i] = sig.fNHighGain[i];
 168       fNLowGain[i]  = sig.fNLowGain[i];
 169     }
 170     fNRef[i] = sig.fNRef[i];
 171   }
 172
 173
 174
 175 }
 176
 177 // assignment operator; use copy ctor to make life easy..
 178 //_____________________________________________________________________
 179 AliCaloCalibSignal& AliCaloCalibSignal::operator = (const AliCaloCalibSignal &source)
 180 {
 181   // assignment operator; use copy ctor
 182   if (&source == this) return *this;
 183
 184   new (this) AliCaloCalibSignal(source);
 185   return *this;
 186 }
 187
 188 //_____________________________________________________________________
 189 void AliCaloCalibSignal::CreateTrees()
 190 {
 191   // initialize trees
 192   // first, regular version
 193   fTreeAmpVsTime = new TTree("fTreeAmpVsTime","Amplitude vs. Time Tree Variables");
 194
 195   fTreeAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I");
 196   fTreeAmpVsTime->Branch("fHour", &fHour, "fHour/D");
 197   fTreeAmpVsTime->Branch("fAmp", &fAmp, "fAmp/D");
 198
 199   // then, average version
 200   fTreeAvgAmpVsTime = new TTree("fTreeAvgAmpVsTime","Average Amplitude vs. Time Tree Variables");
 201
 202   fTreeAvgAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I");
 203   fTreeAvgAmpVsTime->Branch("fHour", &fHour, "fHour/D");
 204   fTreeAvgAmpVsTime->Branch("fAvgAmp", &fAvgAmp, "fAvgAmp/D");
 205   fTreeAvgAmpVsTime->Branch("fRMS", &fRMS, "fRMS/D");
 206
 207   // then same for LED..
 208   fTreeLEDAmpVsTime = new TTree("fTreeLEDAmpVsTime","LED Amplitude vs. Time Tree Variables");
 209   fTreeLEDAmpVsTime->Branch("fRefNum", &fRefNum, "fRefNum/I");
 210   fTreeLEDAmpVsTime->Branch("fHour", &fHour, "fHour/D");
 211   fTreeLEDAmpVsTime->Branch("fAmp", &fAmp, "fAmp/D");
 212
 213   fTreeLEDAvgAmpVsTime = new TTree("fTreeLEDAvgAmpVsTime","Average LED Amplitude vs. Time Tree Variables");
 214   fTreeLEDAvgAmpVsTime->Branch("fRefNum", &fRefNum, "fRefNum/I");
 215   fTreeLEDAvgAmpVsTime->Branch("fHour", &fHour, "fHour/D");
 216   fTreeLEDAvgAmpVsTime->Branch("fAvgAmp", &fAvgAmp, "fAvgAmp/D");
 217   fTreeLEDAvgAmpVsTime->Branch("fRMS", &fRMS, "fRMS/D");
 218
 219   return;
 220 }
 221
 222 //_____________________________________________________________________
 223 void AliCaloCalibSignal::ResetInfo()
 224 { // reset trees and counters
 225   Zero(); // set all counters to 0
 226   DeleteTrees(); // delete previous stuff
 227   CreateTrees(); // and create some new ones
 228   return;
 229 }
 230
 231 //_____________________________________________________________________
 232 void AliCaloCalibSignal::Zero()
 233 {
 234   // set all counters to 0; not cuts etc. though
 235   fHour = 0;
 236   fLatestHour = 0;
 237   fNEvents = 0;
 238   fNAcceptedEvents = 0;
 239
 240   // Set the number of points for each tower: Amp vs. Time
 241   memset(fNHighGain, 0, sizeof(fNHighGain));
 242   memset(fNLowGain, 0, sizeof(fNLowGain));
 243   // and LED reference
 244   memset(fNRef, 0, sizeof(fNRef));
 245
 246   return;
 247 }
 248
 249 //_____________________________________________________________________
 250 Bool_t AliCaloCalibSignal::CheckFractionAboveAmp(const int *iAmpVal,
 251                                                  int resultArray[])
 252 { // check fraction of towers, per column, that are above amplitude cut
 253   Bool_t returnCode = false;
 254
 255   int iTowerNum = 0;
 256   double fraction = 0;
 257   for (int i = 0; i<fModules; i++) {
 258     for (int j = 0; j<fColumns; j++) {
 259       int nAbove = 0;
 260       for (int k = 0; k<fRows; k++) {
 261         iTowerNum = GetTowerNum(i,j,k);
 262         if (iAmpVal[iTowerNum] > fAmpCut) {
 263           nAbove++;
 264         }
 265       }
 266       resultArray[i*fColumns +j] = 0; // init. to denied
 267       if (nAbove > 0) {
 268         fraction = (1.0*nAbove) / fRows;
 269         /*
 270         printf("DS mod %d col %d nAbove %d fraction %3.2f\n",
 271                i, j, nAbove, fraction);
 272         */
 273         if (fraction > fReqFractionAboveAmpCutVal) {
 274           resultArray[i*fColumns + j] = nAbove;
 275           returnCode = true;
 276         }
 277       }
 278     }
 279   } // modules loop
 280
 281   return returnCode;
 282 }
 283
 284
 285 //_____________________________________________________________________
 286 Bool_t AliCaloCalibSignal::CheckLEDRefAboveAmp(const int *iAmpVal,
 287                                                int resultArray[])
 288 { // check which LEDRef/Mon strips are above amplitude cut
 289   Bool_t returnCode = false;
 290
 291   int iRefNum = 0;
 292   int gain = 1; // look at high gain; this should be rather saturated usually..
 293   for (int i = 0; i<fModules; i++) {
 294     for (int j = 0; j<fLEDRefs; j++) {
 295       iRefNum = GetRefNum(i, j, gain);
 296       if (iAmpVal[iRefNum] > fAmpCutLEDRef) {
 297         resultArray[i*fLEDRefs +j] = 1; // enough signal
 298         returnCode = true;
 299       }
 300       else {
 301         resultArray[i*fLEDRefs +j] = 0; // not enough signal
 302       }
 303
 304       /*
 305       printf("DS mod %d LEDRef %d ampVal %d\n",
 306              i, j, iAmpVal[iRefNum]);
 307       */
 308     } // LEDRefs
 309   } // modules loop
 310
 311   return returnCode;
 312 }
 313
 314 // Parameter/cut handling
 315 //_____________________________________________________________________
 316 void AliCaloCalibSignal::SetParametersFromFile(const char *parameterFile)
 317 { // set parameters from file
 318   static const string delimitor("::");
 319
 320   // open, check input file
 321   ifstream in( parameterFile );
 322   if( !in ) {
 323     printf("in AliCaloCalibSignal::SetParametersFromFile - Using default/run_time parameters.\n");
 324     return;
 325   }
 326
 327   // Note: this method is a bit more complicated than it really has to be
 328   // - allowing for multiple entries per line, arbitrary order of the
 329   // different variables etc. But I wanted to try and do this in as
 330   // correct a C++ way as I could (as an exercise).
 331
 332   // read in
 333   char readline[1024];
 334   while ((in.rdstate() & ios::failbit) == 0 ) {
 335
 336     // Read into the raw char array and then construct a string
 337     // to do the searching
 338     in.getline(readline, 1024);
 339     istringstream s(readline);
 340
 341     while ( ( s.rdstate() & ios::failbit ) == 0 ) {
 342
 343       string keyValue;
 344       s >> keyValue;
 345
 346       // check stream status
 347       if( ( s.rdstate() & ios::failbit ) == ios::failbit ) break;
 348
 349       // skip rest of line if comments found
 350       if( keyValue.substr( 0, 2 ) == "//" ) break;
 351
 352       // look for "::" in keyValue pair
 353       size_t position = keyValue.find( delimitor );
 354       if( position == string::npos ) {
 355         printf("wrong format for key::value pair: %s\n", keyValue.c_str());
 356       }
 357
 358       // split keyValue pair
 359       string key( keyValue.substr( 0, position ) );
 360       string value( keyValue.substr( position+delimitor.size(),
 361                                       keyValue.size()-delimitor.size() ) );
 362
 363       // check value does not contain a new delimitor
 364       if( value.find( delimitor ) != string::npos ) {
 365         printf("wrong format for key::value pair: %s\n", keyValue.c_str());
 366       }
 367
 368       // debug: check key value pair
 369       // printf("AliCaloCalibSignal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
 370
 371       // if the key matches with something we expect, we assign the new value
 372       if ( (key == "fAmpCut") || (key == "fReqFractionAboveAmpCutVal") ||
 373            (key == "fAmpCutLEDRef") || (key == "fSecInAverage") ) {
 374         istringstream iss(value);
 375         printf("AliCaloCalibSignal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
 376
 377         if (key == "fAmpCut") {
 378           iss >> fAmpCut;
 379         }
 380         else if (key == "fReqFractionAboveAmpCutVal") {
 381           iss >> fReqFractionAboveAmpCutVal;
 382         }
 383         else if (key == "fAmpCutLEDRef") {
 384           iss >> fAmpCutLEDRef;
 385         }
 386         else if (key == "fSecInAverage") {
 387           iss >> fSecInAverage;
 388         }
 389       } // some match found/expected
 390
 391     }
 392   }
 393
 394   in.close();
 395   return;
 396 }
 397
 398 //_____________________________________________________________________
 399 void AliCaloCalibSignal::WriteParametersToFile(const char *parameterFile)
 400 { // write parameters to file
 401   static const string delimitor("::");
 402   ofstream out( parameterFile );
 403   out << "// " << parameterFile << endl;
 404   out << "fAmpCut" << "::" << fAmpCut << endl;
 405   out << "fReqFractionAboveAmpCutVal" << "::" << fReqFractionAboveAmpCutVal << endl;
 406   out << "fAmpCutLEDRef" << "::" << fAmpCutLEDRef << endl;
 407   out << "fSecInAverage" << "::" << fSecInAverage << endl;
 408
 409   out.close();
 410   return;
 411 }
 412
 413 //_____________________________________________________________________
 414 Bool_t AliCaloCalibSignal::AddInfo(const AliCaloCalibSignal *sig)
 415 {
 416   // note/FIXME: we are not yet adding correctly the info for fN{HighGain,LowGain,Ref} here - but consider this a feature for now (20080905): we'll do Analyze() unless entries were found for a tower in this original object.
 417
 418   // add info from sig's TTrees to ours..
 419   TTree *sigAmp = sig->GetTreeAmpVsTime();
 420   TTree *sigAvgAmp = sig->GetTreeAvgAmpVsTime();
 421
 422   // we could try some merging via TList or what also as a more elegant approach
 423   // but I wanted with the stupid/simple and hopefully safe approach of looping
 424   // over what we want to add..
 425
 426   // associate variables for sigAmp and sigAvgAmp:
 427   sigAmp->SetBranchAddress("fChannelNum",&fChannelNum);
 428   sigAmp->SetBranchAddress("fHour",&fHour);
 429   sigAmp->SetBranchAddress("fAmp",&fAmp);
 430
 431   // loop over the trees.. note that since we use the same variables we should not need
 432   // to do any assignments between the getting and filling
 433   for (int i=0; i<sigAmp->GetEntries(); i++) {
 434     sigAmp->GetEntry(i);
 435     fTreeAmpVsTime->Fill();
 436   }
 437
 438   sigAvgAmp->SetBranchAddress("fChannelNum",&fChannelNum);
 439   sigAvgAmp->SetBranchAddress("fHour",&fHour);
 440   sigAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
 441   sigAvgAmp->SetBranchAddress("fRMS",&fRMS);
 442
 443   for (int i=0; i<sigAvgAmp->GetEntries(); i++) {
 444     sigAvgAmp->GetEntry(i);
 445     fTreeAvgAmpVsTime->Fill();
 446   }
 447
 448   // also LED..
 449   TTree *sigLEDAmp = sig->GetTreeLEDAmpVsTime();
 450   TTree *sigLEDAvgAmp = sig->GetTreeLEDAvgAmpVsTime();
 451
 452   // associate variables for sigAmp and sigAvgAmp:
 453   sigLEDAmp->SetBranchAddress("fRefNum",&fRefNum);
 454   sigLEDAmp->SetBranchAddress("fHour",&fHour);
 455   sigLEDAmp->SetBranchAddress("fAmp",&fAmp);
 456
 457   // loop over the trees.. note that since we use the same variables we should not need
 458   // to do any assignments between the getting and filling
 459   for (int i=0; i<sigLEDAmp->GetEntries(); i++) {
 460     sigLEDAmp->GetEntry(i);
 461     fTreeLEDAmpVsTime->Fill();
 462   }
 463
 464   sigLEDAvgAmp->SetBranchAddress("fRefNum",&fRefNum);
 465   sigLEDAvgAmp->SetBranchAddress("fHour",&fHour);
 466   sigLEDAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
 467   sigLEDAvgAmp->SetBranchAddress("fRMS",&fRMS);
 468
 469   for (int i=0; i<sigLEDAvgAmp->GetEntries(); i++) {
 470     sigLEDAvgAmp->GetEntry(i);
 471     fTreeLEDAvgAmpVsTime->Fill();
 472   }
 473
 474
 475   return kTRUE;//We hopefully succesfully added info from the supplied object
 476 }
 477
 478 //_____________________________________________________________________
 479 Bool_t AliCaloCalibSignal::ProcessEvent(AliRawReader *rawReader)
 480 {
 481   // if fMapping is NULL the rawstream will crate its own mapping
 482   AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
 483   if (fDetType == kEmCal) {
 484     rawReader->Select("EMCAL", 0, AliEMCALGeoParams::fgkLastAltroDDL) ; //select EMCAL DDL range
 485   }
 486   return ProcessEvent( &rawStream, rawReader->GetTimestamp() );
 487 }
 488
 489 //_____________________________________________________________________
 490 Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStreamV3 *in, UInt_t Timestamp)
 491 {
 492   // Method to process=analyze one event in the data stream
 493   if (!in) return kFALSE; //Return right away if there's a null pointer
 494
 495   fNEvents++; // one more event
 496
 497   // use maximum numbers to set array sizes
 498   int iAmpValHighGain[fgkMaxTowers];
 499   int iAmpValLowGain[fgkMaxTowers];
 500   memset(iAmpValHighGain, 0, sizeof(iAmpValHighGain));
 501   memset(iAmpValLowGain, 0, sizeof(iAmpValLowGain));
 502
 503   // also for LED reference
 504   int iLEDAmpVal[fgkMaxRefs * 2]; // factor 2 is for the two gain values
 505   memset(iLEDAmpVal, 0, sizeof(iLEDAmpVal));
 506
 507   int sample = 0; // temporary value
 508   int gain = 0; // high or low gain
 509
 510   // Number of Low and High gain, and LED Ref, channels for this event:
 511   int nLowChan = 0;
 512   int nHighChan = 0;
 513   int nLEDRefChan = 0;
 514
 515   int iTowerNum = 0; // array index for regular towers
 516   int iRefNum = 0; // array index for LED references
 517
 518   // loop first to get the fraction of channels with amplitudes above cut
 519
 520   while (in->NextDDL()) {
 521     while (in->NextChannel()) {
 522
 523       // counters
 524       int max = AliEMCALGeoParams::fgkSampleMin, min = AliEMCALGeoParams::fgkSampleMax; // min and max sample values
 525       int nsamples = 0;
 526
 527       while (in->NextBunch()) {
 528         const UShort_t *sig = in->GetSignals();
 529         nsamples += in->GetBunchLength();
 530         for (Int_t i = 0; i < in->GetBunchLength(); i++) {
 531           sample = sig[i];
 532
 533           // check if it's a min or max value
 534           if (sample < min) min = sample;
 535           if (sample > max) max = sample;
 536
 537         } // loop over samples in bunch
 538       } // loop over bunches
 539
 540       if (nsamples > 0) { // this check is needed for when we have zero-supp. on, but not sparse readout
 541
 542       gain = -1; // init to not valid value
 543       //If we're here then we're done with this tower
 544       if ( in->IsLowGain() ) {
 545         gain = 0;
 546       }
 547       else if ( in->IsHighGain() ) {
 548         gain = 1;
 549       }
 550       else if ( in->IsLEDMonData() ) {
 551         gain = in->GetRow(); // gain coded in (in RCU/Altro mapping) as Row info for LED refs..
 552       }
 553
 554       // it should be enough to check the SuperModule info for each DDL really, but let's keep it here for now
 555       int arrayPos = in->GetModule(); //The modules are numbered starting from 0
 556       //Debug
 557       if (arrayPos < 0 || arrayPos >= fModules) {
 558         printf("AliCaloCalibSignal::ProcessEvent = Oh no: arrayPos = %i.\n", arrayPos);
 559         return kFALSE;
 560       }
 561
 562       if ( in->IsHighGain() || in->IsLowGain() ) { // regular tower
 563         // get tower number for AmpVal array
 564         iTowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow());
 565
 566         if (gain == 0) {
 567           // fill amplitude into the array
 568           iAmpValLowGain[iTowerNum] = max - min;
 569           nLowChan++;
 570         }
 571         else if (gain==1) {//fill the high gain ones
 572           // fill amplitude into the array
 573           iAmpValHighGain[iTowerNum] = max - min;
 574           nHighChan++;
 575         }//end if gain
 576       } // regular tower
 577       else if ( in->IsLEDMonData() ) { // LED ref.;
 578         // strip # is coded is 'column' in the channel maps
 579         iRefNum = GetRefNum(arrayPos, in->GetColumn(), gain);
 580         iLEDAmpVal[iRefNum] = max - min;
 581         nLEDRefChan++;
 582       } // end of LED ref
 583
 584       } // nsamples>0 check, some data found for this channel; not only trailer/header
 585     } // end while over channel
 586
 587   }//end while over DDL's, of input stream
 588
 589   in->Reset(); // just in case the next customer forgets to check if the stream was reset..
 590
 591   // now check if it was an LED event, using the LED Reference info per strip
 592
 593   // by default all columns are accepted (init check to > 0)
 594   int checkResultArray[AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALCols];
 595   for (int ia=0; ia<(AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALCols); ia++) {
 596     checkResultArray[ia] = 1;
 597   }
 598   if (fReqFractionAboveAmp) {
 599     bool ok = false;
 600     if (nHighChan > 0) {
 601       ok = CheckFractionAboveAmp(iAmpValHighGain, checkResultArray);
 602     }
 603     if (!ok) return false; // skip event
 604   }
 605
 606   // by default all columns are accepted (init check to > 0)
 607   int checkResultArrayLEDRef[AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALLEDRefs];
 608   for (int ia=0; ia<(AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALLEDRefs); ia++) {
 609     checkResultArrayLEDRef[ia] = 1;
 610   }
 611   if (fReqLEDRefAboveAmpCutVal) {
 612     bool ok = false;
 613     if (nLEDRefChan > 0) {
 614       ok = CheckLEDRefAboveAmp(iLEDAmpVal, checkResultArrayLEDRef);
 615     }
 616     if (!ok) return false; // skip event
 617   }
 618
 619   fNAcceptedEvents++; // one more event accepted
 620
 621   if (fStartTime == 0) { // if start-timestamp wasn't set,we'll pick it up from the first event we encounter
 622     fStartTime = Timestamp;
 623   }
 624
 625   fHour = (Timestamp - fStartTime)/(double)fgkNumSecInHr;
 626   if (fLatestHour < fHour) {
 627     fLatestHour = fHour;
 628   }
 629
 630   // it is a led event, now fill TTree
 631   // We also do the activity check for LEDRefs/Strips, but need to translate between column
 632   // and strip indices for that; based on these relations:
 633   // iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
 634   // iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2;
 635   // which leads to
 636   // iColFirst = (iSM%2==0) ? iStrip*2 : (AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iStrip)*2;
 637
 638   for(int i=0; i<fModules; i++){
 639     for(int j=0; j<fColumns; j++) {
 640       int iStrip = (i%2==0) ? j/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - j/2;
 641       if (checkResultArray[i*fColumns + j]>0  && checkResultArrayLEDRef[i*fLEDRefs + iStrip]>0) { // column passed check
 642       for(int k=0; k<fRows; k++){
 643
 644         iTowerNum = GetTowerNum(i, j, k);
 645
 646         if(iAmpValHighGain[iTowerNum]) {
 647           fAmp = iAmpValHighGain[iTowerNum];
 648           fChannelNum = GetChannelNum(i,j,k,1);
 649           fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValHighGain[iTowerNum]);
 650           fNHighGain[iTowerNum]++;
 651         }
 652         if(iAmpValLowGain[iTowerNum]) {
 653           fAmp = iAmpValLowGain[iTowerNum];
 654           fChannelNum = GetChannelNum(i,j,k,0);
 655           fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValLowGain[iTowerNum]);
 656           fNLowGain[iTowerNum]++;
 657         }
 658       } // rows
 659       } // column passed check, and LED Ref for strip passed check (if any)
 660     } // columns
 661
 662     // also LED refs
 663     for(int j=0; j<fLEDRefs; j++){
 664       int iColFirst = (i%2==0) ? j*2 : (AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - j)*2; //CHECKME!!!
 665       if ( ((checkResultArray[i*fColumns + iColFirst]>0) || (checkResultArray[i*fColumns + iColFirst + 1]>0)) && // at least one column in strip passed check
 666            (checkResultArrayLEDRef[i*fLEDRefs + j]>0) ) { // and LED Ref passed checks
 667         for (gain=0; gain<2; gain++) {
 668           fRefNum = GetRefNum(i, j, gain);
 669           if (iLEDAmpVal[fRefNum]) {
 670             fAmp = iLEDAmpVal[fRefNum];
 671             fTreeLEDAmpVsTime->Fill();//fRefNum,fHour,fAmp);
 672             fNRef[fRefNum]++;
 673           }
 674         } // gain
 675       } // at least one column in strip passed check, and LED Ref passed check (if any)
 676     }
 677
 678   } // modules
 679
 680   return kTRUE;
 681 }
 682
 683 //_____________________________________________________________________
 684 Bool_t AliCaloCalibSignal::Save(TString fileName)
 685 {
 686   //Saves all the TTrees to the designated file
 687
 688   TFile destFile(fileName, "recreate");
 689
 690   if (destFile.IsZombie()) {
 691     return kFALSE;
 692   }
 693
 694   destFile.cd();
 695
 696   // save the trees
 697   fTreeAmpVsTime->Write();
 698   fTreeLEDAmpVsTime->Write();
 699   if (fUseAverage) {
 700     Analyze(); // get the latest and greatest averages
 701     fTreeAvgAmpVsTime->Write();
 702     fTreeLEDAvgAmpVsTime->Write();
 703   }
 704
 705   destFile.Close();
 706
 707   return kTRUE;
 708 }
 709
 710 //_____________________________________________________________________
 711 Bool_t AliCaloCalibSignal::Analyze()
 712 {
 713   // Fill the tree holding the average values
 714   if (!fUseAverage) { return kFALSE; }
 715
 716   // Reset the average TTree if Analyze has already been called earlier,
 717   // meaning that the TTree could have been partially filled
 718   if (fTreeAvgAmpVsTime->GetEntries() > 0) {
 719     fTreeAvgAmpVsTime->Reset();
 720   }
 721
 722   //0: setup variables for the TProfile plots that we'll use to do the averages
 723   int numProfBins = 0;
 724   double timeMin = 0;
 725   double timeMax = 0;
 726   if (fSecInAverage > 0) {
 727     numProfBins = (int)( (fLatestHour*fgkNumSecInHr)/fSecInAverage + 1 ); // round-off
 728   }
 729   numProfBins += 2; // add extra buffer : first and last
 730   double binSize = 1.0*fSecInAverage / fgkNumSecInHr;
 731   timeMin = - binSize;
 732   timeMax = timeMin + numProfBins*binSize;
 733
 734   //1: set up TProfiles for the towers that had data
 735   TProfile * profile[fgkMaxTowers*2]; // *2 is since we include both high and low gains
 736   memset(profile, 0, sizeof(profile));
 737   const Int_t buffersize = 200;
 738   char name[buffersize]; // for profile id and title
 739   int iTowerNum = 0;
 740
 741   for (int i = 0; i<fModules; i++) {
 742     for (int ic=0; ic<fColumns; ic++){
 743       for (int ir=0; ir<fRows; ir++) {
 744
 745         iTowerNum = GetTowerNum(i, ic, ir);
 746         // high gain
 747         if (fNHighGain[iTowerNum] > 0) {
 748           fChannelNum = GetChannelNum(i, ic, ir, 1);
 749           snprintf(name,buffersize,"profileChan%d", fChannelNum);
 750           profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
 751         }
 752
 753         // same for low gain
 754         if (fNLowGain[iTowerNum] > 0) {
 755           fChannelNum = GetChannelNum(i, ic, ir, 0);
 756           snprintf(name,buffersize,"profileChan%d", fChannelNum);
 757           profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
 758         }
 759
 760       } // rows
 761     } // columns
 762   } // modules
 763
 764   //2: fill profiles by looping over tree
 765   // Set addresses for tree-readback also
 766   fTreeAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
 767   fTreeAmpVsTime->SetBranchAddress("fHour", &fHour);
 768   fTreeAmpVsTime->SetBranchAddress("fAmp", &fAmp);
 769
 770   for (int ient=0; ient<fTreeAmpVsTime->GetEntries(); ient++) {
 771     fTreeAmpVsTime->GetEntry(ient);
 772     if (profile[fChannelNum]) {
 773       // profile should always have been created above, for active channels
 774       profile[fChannelNum]->Fill(fHour, fAmp);
 775     }
 776   }
 777
 778   // re-associating the branch addresses here seems to be needed for OK 'average' storage
 779   fTreeAvgAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
 780   fTreeAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
 781   fTreeAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
 782   fTreeAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);
 783
 784   //3: fill avg tree by looping over the profiles
 785   for (fChannelNum = 0; fChannelNum<(fgkMaxTowers*2); fChannelNum++) {
 786     if (profile[fChannelNum]) { // profile was created
 787       if (profile[fChannelNum]->GetEntries() > 0) { // profile had some entries
 788         for(int it=0; it<numProfBins; it++) {
 789           if (profile[fChannelNum]->GetBinEntries(it+1) > 0) {
 790             fAvgAmp = profile[fChannelNum]->GetBinContent(it+1);
 791             fHour = profile[fChannelNum]->GetBinCenter(it+1);
 792             fRMS = profile[fChannelNum]->GetBinError(it+1);
 793             fTreeAvgAmpVsTime->Fill();
 794           } // some entries for this bin
 795         } // loop over bins
 796       } // some entries for this profile
 797     } // profile exists
 798   } // loop over all possible channels
 799
 800
 801   // and finally, go through same exercise for LED also..
 802
 803   //1: set up TProfiles for the towers that had data
 804   TProfile * profileLED[fgkMaxRefs*2]; // *2 is since we include both high and low gains
 805   memset(profileLED, 0, sizeof(profileLED));
 806
 807   for (int i = 0; i<fModules; i++) {
 808     for(int j=0; j<fLEDRefs; j++){
 809       for (int gain=0; gain<2; gain++) {
 810         fRefNum = GetRefNum(i, j, gain);
 811         if (fNRef[fRefNum] > 0) {
 812           sprintf(name, "profileLEDRef%d", fRefNum);
 813           profileLED[fRefNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
 814         }
 815       }// gain
 816     }
 817   } // modules
 818
 819   //2: fill profiles by looping over tree
 820   // Set addresses for tree-readback also
 821   fTreeLEDAmpVsTime->SetBranchAddress("fRefNum", &fRefNum);
 822   fTreeLEDAmpVsTime->SetBranchAddress("fHour", &fHour);
 823   fTreeLEDAmpVsTime->SetBranchAddress("fAmp", &fAmp);
 824
 825   for (int ient=0; ient<fTreeLEDAmpVsTime->GetEntries(); ient++) {
 826     fTreeLEDAmpVsTime->GetEntry(ient);
 827     if (profileLED[fRefNum]) {
 828       // profile should always have been created above, for active channels
 829       profileLED[fRefNum]->Fill(fHour, fAmp);
 830     }
 831   }
 832
 833   // re-associating the branch addresses here seems to be needed for OK 'average' storage
 834   fTreeLEDAvgAmpVsTime->SetBranchAddress("fRefNum", &fRefNum);
 835   fTreeLEDAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
 836   fTreeLEDAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
 837   fTreeLEDAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);
 838
 839   //3: fill avg tree by looping over the profiles
 840   for (fRefNum = 0; fRefNum<(fgkMaxRefs*2); fRefNum++) {
 841     if (profileLED[fRefNum]) { // profile was created
 842       if (profileLED[fRefNum]->GetEntries() > 0) { // profile had some entries
 843         for(int it=0; it<numProfBins; it++) {
 844           if (profileLED[fRefNum]->GetBinEntries(it+1) > 0) {
 845             fAvgAmp = profileLED[fRefNum]->GetBinContent(it+1);
 846             fHour = profileLED[fRefNum]->GetBinCenter(it+1);
 847             fRMS = profileLED[fRefNum]->GetBinError(it+1);
 848             fTreeLEDAvgAmpVsTime->Fill();
 849           } // some entries for this bin
 850         } // loop over bins
 851       } // some entries for this profile
 852     } // profile exists
 853   } // loop over all possible channels
 854
 855   // OK, we're done..
 856
 857   return kTRUE;
 858 }
 859
 860 //_____________________________________________________________________
 861 Bool_t AliCaloCalibSignal::DecodeChannelNum(const int chanId,
 862                                             int *imod, int *icol, int *irow, int *igain) const
 863 { // return the module, column, row, and gain for a given channel number
 864   *igain = chanId/(fModules*fColumns*fRows);
 865   *imod = (chanId/(fColumns*fRows)) % fModules;
 866   *icol = (chanId/fRows) % fColumns;
 867   *irow = chanId % fRows;
 868   return kTRUE;
 869 }
 870
 871 //_____________________________________________________________________
 872 Bool_t AliCaloCalibSignal::DecodeRefNum(const int refId,
 873                                         int *imod, int *istripMod, int *igain) const
 874 { // return the module, stripModule, and gain for a given reference number
 875   *igain = refId/(fModules*fLEDRefs);
 876   *imod = (refId/(fLEDRefs)) % fModules;
 877   *istripMod = refId % fLEDRefs;
 878   return kTRUE;
 879 }