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