1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 $ */
17 //________________________________________________________________________
19 // A help class for monitoring and calibration tools: MOOD, AMORE etc.,
20 // It can be created and used a la (ctor):
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();
28 // fSignals->ProcessEvent(fCaloRawStream,fRawEventHeaderBase);
30 // fSignals->GetXXX..()
32 //________________________________________________________________________
40 #include "AliRawReader.h"
41 #include "AliCaloRawStreamV3.h"
43 #include "AliCaloConstants.h"
44 #include "AliCaloBunchInfo.h"
45 #include "AliCaloFitResults.h"
46 #include "AliCaloRawAnalyzer.h"
47 #include "AliCaloRawAnalyzerFactory.h"
50 #include "AliCaloCalibSignal.h"
52 ClassImp(AliCaloCalibSignal)
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;
63 // ctor; initialize everything in order to avoid compiler warnings
64 // put some reasonable defaults
65 AliCaloCalibSignal::AliCaloCalibSignal(kDetType detectorType) :
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),
91 fTreeAvgAmpVsTime(NULL),
92 fTreeLEDAmpVsTime(NULL),
93 fTreeLEDAvgAmpVsTime(NULL)
95 //Default constructor. First we set the detector-type related constants.
96 if (detectorType == kPhos) {
97 fColumns = fgkPhosCols;
99 fLEDRefs = fgkPhosLEDRefs;
100 fModules = fgkPhosModules;
101 fCaloString = "PHOS";
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";
114 fDetType = detectorType;
115 SetFittingAlgorithm(Algo::kStandard);
116 ResetInfo(); // trees and counters
120 //_____________________________________________________________________
121 AliCaloCalibSignal::~AliCaloCalibSignal()
126 //_____________________________________________________________________
127 void AliCaloCalibSignal::DeleteTrees()
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;
144 //_____________________________________________________________________
145 //AliCaloCalibSignal::AliCaloCalibSignal(const AliCaloCalibSignal &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()),
169 // fTreeAvgAmpVsTime(),
170 // fTreeLEDAmpVsTime(),
171 // fTreeLEDAvgAmpVsTime()
173 // // also the TTree contents
175 // for (Int_t i = 0; i<fgkMaxTowers; i++) {
176 // fNHighGain[i] = sig.fNHighGain[i];
177 // fNLowGain[i] = sig.fNLowGain[i];
179 // for (Int_t i = 0; i<(2*fgkMaxRefs); i++) {
180 // fNRef[i] = sig.fNRef[i];
186 // assignment operator; use copy ctor to make life easy..
187 //_____________________________________________________________________
188 //AliCaloCalibSignal& AliCaloCalibSignal::operator = (const AliCaloCalibSignal &source)
190 // // assignment operator; use copy ctor
191 // if (&source == this) return *this;
193 // new (this) AliCaloCalibSignal(source);
197 //_____________________________________________________________________
198 void AliCaloCalibSignal::CreateTrees()
201 // first, regular version
202 fTreeAmpVsTime = new TTree("fTreeAmpVsTime","Amplitude vs. Time Tree Variables");
204 fTreeAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I");
205 fTreeAmpVsTime->Branch("fHour", &fHour, "fHour/D");
206 fTreeAmpVsTime->Branch("fAmp", &fAmp, "fAmp/D");
208 // then, average version
209 fTreeAvgAmpVsTime = new TTree("fTreeAvgAmpVsTime","Average Amplitude vs. Time Tree Variables");
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");
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");
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");
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
240 //_____________________________________________________________________
241 void AliCaloCalibSignal::Zero()
243 // set all counters to 0; not cuts etc. though
247 fNAcceptedEvents = 0;
249 // Set the number of points for each tower: Amp vs. Time
250 memset(fNHighGain, 0, sizeof(fNHighGain));
251 memset(fNLowGain, 0, sizeof(fNLowGain));
253 memset(fNRef, 0, sizeof(fNRef));
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;
266 for (int i = 0; i<fModules; i++) {
267 for (int j = 0; j<fColumns; j++) {
269 for (int k = 0; k<fRows; k++) {
270 iTowerNum = GetTowerNum(i,j,k);
271 if (iAmpVal[iTowerNum] > fAmpCut) {
275 resultArray[i*fColumns +j] = 0; // init. to denied
277 fraction = (1.0*nAbove) / fRows;
279 printf("DS mod %d col %d nAbove %d fraction %3.2f\n",
280 i, j, nAbove, fraction);
282 if (fraction > fReqFractionAboveAmpCutVal) {
283 resultArray[i*fColumns + j] = nAbove;
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;
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
310 resultArray[i*fLEDRefs +j] = 0; // not enough signal
314 printf("DS mod %d LEDRef %d ampVal %d\n",
315 i, j, iAmpVal[iRefNum]);
323 // Parameter/cut handling
324 //_____________________________________________________________________
325 void AliCaloCalibSignal::SetParametersFromFile(const char *parameterFile)
326 { // set parameters from file
327 static const string delimitor("::");
329 // open, check input file
330 ifstream in( parameterFile );
332 printf("in AliCaloCalibSignal::SetParametersFromFile - Using default/run_time parameters.\n");
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).
343 while ((in.rdstate() & ios::failbit) == 0 ) {
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);
350 while ( ( s.rdstate() & ios::failbit ) == 0 ) {
355 // check stream status
356 if( ( s.rdstate() & ios::failbit ) == ios::failbit ) break;
358 // skip rest of line if comments found
359 if( keyValue.substr( 0, 2 ) == "//" ) break;
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());
367 // split keyValue pair
368 string key( keyValue.substr( 0, position ) );
369 string value( keyValue.substr( position+delimitor.size(),
370 keyValue.size()-delimitor.size() ) );
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());
377 // debug: check key value pair
378 // printf("AliCaloCalibSignal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
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());
387 if (key == "fAmpCut") {
390 else if (key == "fReqFractionAboveAmpCutVal") {
391 iss >> fReqFractionAboveAmpCutVal;
393 else if (key == "fAmpCutLEDRef") {
394 iss >> fAmpCutLEDRef;
396 else if (key == "fSecInAverage") {
397 iss >> fSecInAverage;
399 else if (key == "fFittingAlgorithm") {
400 iss >> fFittingAlgorithm;
401 SetFittingAlgorithm( fFittingAlgorithm );
403 else if (key == "fDownscale") {
406 } // some match found/expected
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;
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
441 //_____________________________________________________________________
442 Bool_t AliCaloCalibSignal::AddInfo(const AliCaloCalibSignal *sig)
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.
446 // add info from sig's TTrees to ours..
447 TTree *sigAmp = sig->GetTreeAmpVsTime();
448 TTree *sigAvgAmp = sig->GetTreeAvgAmpVsTime();
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..
454 // associate variables for sigAmp and sigAvgAmp:
455 sigAmp->SetBranchAddress("fChannelNum",&fChannelNum);
456 sigAmp->SetBranchAddress("fHour",&fHour);
457 sigAmp->SetBranchAddress("fAmp",&fAmp);
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++) {
463 fTreeAmpVsTime->Fill();
466 sigAvgAmp->SetBranchAddress("fChannelNum",&fChannelNum);
467 sigAvgAmp->SetBranchAddress("fHour",&fHour);
468 sigAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
469 sigAvgAmp->SetBranchAddress("fRMS",&fRMS);
471 for (int i=0; i<sigAvgAmp->GetEntries(); i++) {
472 sigAvgAmp->GetEntry(i);
473 fTreeAvgAmpVsTime->Fill();
477 TTree *sigLEDAmp = sig->GetTreeLEDAmpVsTime();
478 TTree *sigLEDAvgAmp = sig->GetTreeLEDAvgAmpVsTime();
480 // associate variables for sigAmp and sigAvgAmp:
481 sigLEDAmp->SetBranchAddress("fRefNum",&fRefNum);
482 sigLEDAmp->SetBranchAddress("fHour",&fHour);
483 sigLEDAmp->SetBranchAddress("fAmp",&fAmp);
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();
492 sigLEDAvgAmp->SetBranchAddress("fRefNum",&fRefNum);
493 sigLEDAvgAmp->SetBranchAddress("fHour",&fHour);
494 sigLEDAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
495 sigLEDAvgAmp->SetBranchAddress("fRMS",&fRMS);
497 for (int i=0; i<sigLEDAvgAmp->GetEntries(); i++) {
498 sigLEDAvgAmp->GetEntry(i);
499 fTreeLEDAvgAmpVsTime->Fill();
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();
521 return kTRUE;//We hopefully succesfully added info from the supplied object
524 //_____________________________________________________________________
525 Bool_t AliCaloCalibSignal::ProcessEvent(AliRawReader *rawReader)
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
533 return ProcessEvent( &rawStream, rawReader->GetTimestamp() );
536 //_____________________________________________________________________
537 Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStreamV3 *in, UInt_t Timestamp)
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
542 fNEvents++; // one more event
544 if ( (fNEvents%fDownscale)!=0 ) return kFALSE; // mechanism to skip some of the input events, if we want
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));
552 // also for LED reference
553 int iLEDAmpVal[fgkMaxRefs * 2]; // factor 2 is for the two gain values
554 memset(iLEDAmpVal, 0, sizeof(iLEDAmpVal));
556 int gain = 0; // high or low gain
558 // Number of Low and High gain, and LED Ref, channels for this event:
563 int iTowerNum = 0; // array index for regular towers
564 int iRefNum = 0; // array index for LED references
566 // loop first to get the fraction of channels with amplitudes above cut
568 while (in->NextDDL()) {
569 while (in->NextChannel()) {
571 vector<AliCaloBunchInfo> bunchlist;
572 while (in->NextBunch()) {
573 bunchlist.push_back( AliCaloBunchInfo(in->GetStartTimeBin(), in->GetBunchLength(), in->GetSignals() ) );
575 if (bunchlist.size() == 0) continue;
577 gain = -1; // init to not valid value
578 //If we're here then we're done with this tower
579 if ( in->IsLowGain() ) {
582 else if ( in->IsHighGain() ) {
585 else if ( in->IsLEDMonData() ) {
586 gain = in->GetRow(); // gain coded in (in RCU/Altro mapping) as Row info for LED refs..
588 else { continue; } // don't try to fit TRU..
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
593 if (arrayPos < 0 || arrayPos >= fModules) {
594 printf("AliCaloCalibSignal::ProcessEvent = Oh no: arrayPos = %i.\n", arrayPos);
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());
604 // fill amplitude into the array
605 iAmpValLowGain[iTowerNum] = (int) res.GetAmp();
608 else if (gain==1) {//fill the high gain ones
609 // fill amplitude into the array
610 iAmpValHighGain[iTowerNum] = (int) res.GetAmp();
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();
621 } // end while over channel
623 }//end while over DDL's, of input stream
625 in->Reset(); // just in case the next customer forgets to check if the stream was reset..
627 // now check if it was an LED event, using the LED Reference info per strip
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;
634 if (fReqFractionAboveAmp) {
637 ok = CheckFractionAboveAmp(iAmpValHighGain, checkResultArray);
639 if (!ok) return false; // skip event
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;
647 if (fReqLEDRefAboveAmpCutVal) {
649 if (nLEDRefChan > 0) {
650 ok = CheckLEDRefAboveAmp(iLEDAmpVal, checkResultArrayLEDRef);
652 if (!ok) return false; // skip event
655 fNAcceptedEvents++; // one more event accepted
657 if (fStartTime == 0) { // if start-timestamp wasn't set,we'll pick it up from the first event we encounter
658 fStartTime = Timestamp;
661 fHour = (Timestamp - fStartTime)/(double)fgkNumSecInHr;
662 if (fLatestHour < fHour) {
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;
672 // iColFirst = (iSM%2==0) ? iStrip*2 : (AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iStrip)*2;
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++){
680 iTowerNum = GetTowerNum(i, j, k);
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]++;
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]++;
695 } // column passed check, and LED Ref for strip passed check (if any)
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);
711 } // at least one column in strip passed check, and LED Ref passed check (if any)
719 //_____________________________________________________________________
720 Bool_t AliCaloCalibSignal::Save(TString fileName)
722 //Saves all the TTrees to the designated file
724 TFile destFile(fileName, "recreate");
726 if (destFile.IsZombie()) {
733 fTreeAmpVsTime->Write();
734 fTreeLEDAmpVsTime->Write();
736 Analyze(); // get the latest and greatest averages
737 fTreeAvgAmpVsTime->Write();
738 fTreeLEDAvgAmpVsTime->Write();
746 //_____________________________________________________________________
747 Bool_t AliCaloCalibSignal::Analyze()
749 // Fill the tree holding the average values
750 if (!fUseAverage) { return kFALSE; }
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();
758 //0: setup variables for the TProfile plots that we'll use to do the averages
762 if (fSecInAverage > 0) {
763 numProfBins = (int)( (fLatestHour*fgkNumSecInHr)/fSecInAverage + 1 ); // round-off
765 numProfBins += 2; // add extra buffer : first and last
766 double binSize = 1.0*fSecInAverage / fgkNumSecInHr;
768 timeMax = timeMin + numProfBins*binSize;
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
777 for (int i = 0; i<fModules; i++) {
778 for (int ic=0; ic<fColumns; ic++){
779 for (int ir=0; ir<fRows; ir++) {
781 iTowerNum = GetTowerNum(i, ic, ir);
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");
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");
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);
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);
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);
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
832 } // some entries for this profile
834 } // loop over all possible channels
837 // and finally, go through same exercise for LED also..
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));
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");
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);
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);
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);
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
887 } // some entries for this profile
889 } // loop over all possible channels
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;
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;