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 **************************************************************************/
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 fPedestals = new AliCaloCalibPedestal( fDetType );
24 // AliCaloCalibPedestal knows how many modules we have for PHOS or EMCAL
25 fNumModules = fPedestals->GetModules();
28 // fPedestals->ProcessEvent(fCaloRawStream);
29 // asked to draw histograms:
30 // fPedestals->GetDeadMap(i)->Draw("col");
32 // fPedestals->GetPeakProfileHighGainRatio((i < fNumModules) ? i : fVisibleModule)->Draw("colz");
34 // The pseudo-code examples above were from the first implementation in MOOD (summer 2007).
35 //________________________________________________________________________
37 //#include "TCanvas.h"
47 #include "AliRawReader.h"
48 #include "AliCaloRawStreamV3.h"
51 #include "AliCaloCalibPedestal.h"
53 ClassImp(AliCaloCalibPedestal)
57 // ctor; initialize everything in order to avoid compiler warnings
58 AliCaloCalibPedestal::AliCaloCalibPedestal(kDetType detectorType) :
62 fPedestalLEDRefLowGain(),
63 fPedestalLEDRefHighGain(),
64 fPeakMinusPedLowGain(),
65 fPeakMinusPedHighGain(),
66 fPeakMinusPedHighGainHisto(),
67 fPedestalLowGainDiff(),
68 fPedestalHighGainDiff(),
69 fPedestalLEDRefLowGainDiff(),
70 fPedestalLEDRefHighGainDiff(),
71 fPeakMinusPedLowGainDiff(),
72 fPeakMinusPedHighGainDiff(),
73 fPedestalLowGainRatio(),
74 fPedestalHighGainRatio(),
75 fPedestalLEDRefLowGainRatio(),
76 fPedestalLEDRefHighGainRatio(),
77 fPeakMinusPedLowGainRatio(),
78 fPeakMinusPedHighGainRatio(),
84 fResurrectedTowers(0),
97 fSelectPedestalSamples(kTRUE),
98 fFirstPedestalSample(0),
99 fLastPedestalSample(15),
101 fWarningThreshold(50),
102 fWarningFraction(0.002),
105 //Default constructor. First we set the detector-type related constants.
106 if (detectorType == kPhos) {
107 fColumns = fgkPhosCols;
109 fLEDRefs = fgkPhosLEDRefs;
110 fModules = fgkPhosModules;
111 fCaloString = "PHOS";
117 //We'll just trust the enum to keep everything in line, so that if detectorType
118 //isn't kPhos then it is kEmCal. Note, however, that this is not necessarily the
119 //case, if someone intentionally gives another number
120 fColumns = AliEMCALGeoParams::fgkEMCALCols;
121 fRows = AliEMCALGeoParams::fgkEMCALRows;
122 fLEDRefs = AliEMCALGeoParams::fgkEMCALLEDRefs;
123 fModules = AliEMCALGeoParams::fgkEMCALModules;
124 fCaloString = "EMCAL";
129 fDetType = detectorType;
131 //Then, loop for the requested number of modules
133 for (int i = 0; i < fModules; i++) {
134 //Pedestals, low gain
135 name = "hPedlowgain";
137 title = "Pedestals, low gain, module ";
139 fPedestalLowGain.Add(new TProfile2D(name, title,
140 fColumns, 0.0, fColumns,
141 fRows, fRowMin, fRowMax,"s"));
143 //Pedestals, high gain
144 name = "hPedhighgain";
146 title = "Pedestals, high gain, module ";
148 fPedestalHighGain.Add(new TProfile2D(name, title,
149 fColumns, 0.0, fColumns,
150 fRows, fRowMin, fRowMax,"s"));
152 //LED Ref/Mon pedestals, low gain
153 name = "hPedestalLEDReflowgain";
155 title = "Pedestal LEDRef, low gain, module ";
157 fPedestalLEDRefLowGain.Add(new TProfile(name, title,
158 fLEDRefs, 0.0, fLEDRefs, "s"));
160 //LED Ref/Mon pedestals, high gain
161 name = "hPedestalLEDRefhighgain";
163 title = "Pedestal LEDRef, high gain, module ";
165 fPedestalLEDRefHighGain.Add(new TProfile(name, title,
166 fLEDRefs, 0.0, fLEDRefs, "s"));
168 //Peak-Pedestals, low gain
169 name = "hPeakMinusPedlowgain";
171 title = "Peak-Pedestal, low gain, module ";
173 fPeakMinusPedLowGain.Add(new TProfile2D(name, title,
174 fColumns, 0.0, fColumns,
175 fRows, fRowMin, fRowMax,"s"));
177 //Peak-Pedestals, high gain
178 name = "hPeakMinusPedhighgain";
180 title = "Peak-Pedestal, high gain, module ";
182 fPeakMinusPedHighGain.Add(new TProfile2D(name, title,
183 fColumns, 0.0, fColumns,
184 fRows, fRowMin, fRowMax,"s"));
186 //Peak-Pedestals, high gain - TH2F histo
187 name = "hPeakMinusPedhighgainHisto";
189 title = "Peak-Pedestal, high gain, module ";
191 fPeakMinusPedHighGainHisto.Add(new TH2F(name, title,
192 fColumns*fRows, 0.0, fColumns*fRows,
197 title = "Dead map, module ";
199 fDeadMap.Add(new TH2D(name, title, fColumns, 0.0, fColumns,
200 fRows, fRowMin, fRowMax));
202 }//end for nModules create the histograms
204 //Compress the arrays, in order to remove the empty objects (a 16 slot array is created by default)
205 fPedestalLowGain.Compress();
206 fPedestalHighGain.Compress();
207 fPedestalLEDRefLowGain.Compress();
208 fPedestalLEDRefHighGain.Compress();
209 fPeakMinusPedLowGain.Compress();
210 fPeakMinusPedHighGain.Compress();
211 fPeakMinusPedHighGainHisto.Compress();
214 // set owner ship for everyone
215 fPedestalLowGain.SetOwner(kTRUE);
216 fPedestalHighGain.SetOwner(kTRUE);
217 fPedestalLEDRefLowGain.SetOwner(kTRUE);
218 fPedestalLEDRefHighGain.SetOwner(kTRUE);
219 fPeakMinusPedLowGain.SetOwner(kTRUE);
220 fPeakMinusPedHighGain.SetOwner(kTRUE);
221 fPeakMinusPedHighGainHisto.SetOwner(kTRUE);
222 fPedestalLowGainDiff.SetOwner(kTRUE);
223 fPedestalHighGainDiff.SetOwner(kTRUE);
224 fPedestalLEDRefLowGainDiff.SetOwner(kTRUE);
225 fPedestalLEDRefHighGainDiff.SetOwner(kTRUE);
226 fPeakMinusPedLowGainDiff.SetOwner(kTRUE);
227 fPeakMinusPedHighGainDiff.SetOwner(kTRUE);
228 fPedestalLowGainRatio.SetOwner(kTRUE);
229 fPedestalHighGainRatio.SetOwner(kTRUE);
230 fPedestalLEDRefLowGainRatio.SetOwner(kTRUE);
231 fPedestalLEDRefHighGainRatio.SetOwner(kTRUE);
232 fPeakMinusPedLowGainRatio.SetOwner(kTRUE);
233 fPeakMinusPedHighGainRatio.SetOwner(kTRUE);
234 fDeadMap.SetOwner(kTRUE);
239 //_____________________________________________________________________
240 AliCaloCalibPedestal::~AliCaloCalibPedestal()
242 if (fReference) delete fReference;//Delete the reference object, if it has been loaded
244 // delete also TObjArray's
245 fPedestalLowGain.Delete();
246 fPedestalHighGain.Delete();
247 fPedestalLEDRefLowGain.Delete();
248 fPedestalLEDRefHighGain.Delete();
249 fPeakMinusPedLowGain.Delete();
250 fPeakMinusPedHighGain.Delete();
251 fPeakMinusPedHighGainHisto.Delete();
252 fPedestalLowGainDiff.Delete();
253 fPedestalHighGainDiff.Delete();
254 fPedestalLEDRefLowGainDiff.Delete();
255 fPedestalLEDRefHighGainDiff.Delete();
256 fPeakMinusPedLowGainDiff.Delete();
257 fPeakMinusPedHighGainDiff.Delete();
258 fPedestalLowGainRatio.Delete();
259 fPedestalHighGainRatio.Delete();
260 fPedestalLEDRefLowGainRatio.Delete();
261 fPedestalLEDRefHighGainRatio.Delete();
262 fPeakMinusPedLowGainRatio.Delete();
263 fPeakMinusPedHighGainRatio.Delete();
269 //_____________________________________________________________________
270 AliCaloCalibPedestal::AliCaloCalibPedestal(const AliCaloCalibPedestal &ped) :
274 fPedestalLEDRefLowGain(),
275 fPedestalLEDRefHighGain(),
276 fPeakMinusPedLowGain(),
277 fPeakMinusPedHighGain(),
278 fPeakMinusPedHighGainHisto(),
279 fPedestalLowGainDiff(),
280 fPedestalHighGainDiff(),
281 fPedestalLEDRefLowGainDiff(),
282 fPedestalLEDRefHighGainDiff(),
283 fPeakMinusPedLowGainDiff(),
284 fPeakMinusPedHighGainDiff(),
285 fPedestalLowGainRatio(),
286 fPedestalHighGainRatio(),
287 fPedestalLEDRefLowGainRatio(),
288 fPedestalLEDRefHighGainRatio(),
289 fPeakMinusPedLowGainRatio(),
290 fPeakMinusPedHighGainRatio(),
292 fNEvents(ped.GetNEvents()),
293 fNChanFills(ped.GetNChanFills()),
294 fDeadTowers(ped.GetDeadTowerCount()),
295 fNewDeadTowers(ped.GetDeadTowerNew()),
296 fResurrectedTowers(ped.GetDeadTowerResurrected()),
297 fReference( 0 ), //! note that we do not try to copy the reference info here
298 fDetType(ped.GetDetectorType()),
299 fColumns(ped.GetColumns()),
300 fRows(ped.GetRows()),
301 fLEDRefs(ped.GetLEDRefs()),
302 fModules(ped.GetModules()),
303 fRowMin(ped.GetRowMin()),
304 fRowMax(ped.GetRowMax()),
305 fRowMultiplier(ped.GetRowMultiplier()),
306 fCaloString(ped.GetCaloString()),
307 fMapping(NULL), //! note that we are not copying the map info
308 fRunNumber(ped.GetRunNumber()),
309 fSelectPedestalSamples(ped.GetSelectPedestalSamples()),
310 fFirstPedestalSample(ped.GetFirstPedestalSample()),
311 fLastPedestalSample(ped.GetLastPedestalSample()),
312 fDeadThreshold(ped.GetDeadThreshold()),
313 fWarningThreshold(ped.GetWarningThreshold()),
314 fWarningFraction(ped.GetWarningFraction()),
315 fHotSigma(ped.GetHotSigma())
317 // Then the ObjArray ones; we add the histograms rather than trying TObjArray = assignment
318 //DS: this has not really been tested yet..
319 for (int i = 0; i < fModules; i++) {
320 fPedestalLowGain.Add( ped.GetPedProfileLowGain(i) );
321 fPedestalHighGain.Add( ped.GetPedProfileHighGain(i) );
322 fPedestalLEDRefLowGain.Add( ped.GetPedLEDRefProfileLowGain(i) );
323 fPedestalLEDRefHighGain.Add( ped.GetPedLEDRefProfileHighGain(i) );
324 fPeakMinusPedLowGain.Add( ped.GetPeakProfileLowGain(i) );
325 fPeakMinusPedHighGain.Add( ped.GetPeakProfileHighGain(i) );
326 fPeakMinusPedHighGainHisto.Add( ped.GetPeakHighGainHisto(i) );
328 fDeadMap.Add( ped.GetDeadMap(i) );
331 //Compress the arrays, in order to remove the empty objects (a 16 slot array is created by default)
332 fPedestalLowGain.Compress();
333 fPedestalHighGain.Compress();
334 fPedestalLEDRefLowGain.Compress();
335 fPedestalLEDRefHighGain.Compress();
336 fPeakMinusPedLowGain.Compress();
337 fPeakMinusPedHighGain.Compress();
338 fPeakMinusPedHighGainHisto.Compress();
341 // set owner ship for everyone
342 fPedestalLowGain.SetOwner(kTRUE);
343 fPedestalHighGain.SetOwner(kTRUE);
344 fPedestalLEDRefLowGain.SetOwner(kTRUE);
345 fPedestalLEDRefHighGain.SetOwner(kTRUE);
346 fPeakMinusPedLowGain.SetOwner(kTRUE);
347 fPeakMinusPedHighGain.SetOwner(kTRUE);
348 fPeakMinusPedHighGainHisto.SetOwner(kTRUE);
349 fPedestalLowGainDiff.SetOwner(kTRUE);
350 fPedestalHighGainDiff.SetOwner(kTRUE);
351 fPedestalLEDRefLowGainDiff.SetOwner(kTRUE);
352 fPedestalLEDRefHighGainDiff.SetOwner(kTRUE);
353 fPeakMinusPedLowGainDiff.SetOwner(kTRUE);
354 fPeakMinusPedHighGainDiff.SetOwner(kTRUE);
355 fPedestalLowGainRatio.SetOwner(kTRUE);
356 fPedestalHighGainRatio.SetOwner(kTRUE);
357 fPedestalLEDRefLowGainRatio.SetOwner(kTRUE);
358 fPedestalLEDRefHighGainRatio.SetOwner(kTRUE);
359 fPeakMinusPedLowGainRatio.SetOwner(kTRUE);
360 fPeakMinusPedHighGainRatio.SetOwner(kTRUE);
361 fDeadMap.SetOwner(kTRUE);
365 // assignment operator; use copy ctor to make life easy..
366 //_____________________________________________________________________
367 AliCaloCalibPedestal& AliCaloCalibPedestal::operator = (const AliCaloCalibPedestal &source)
369 // assignment operator; use copy ctor
370 if (&source == this) return *this;
372 new (this) AliCaloCalibPedestal(source);
376 //_____________________________________________________________________
377 void AliCaloCalibPedestal::Reset()
379 // Reset all arrays/histograms
380 for (int i = 0; i < fModules; i++) {
381 GetPedProfileLowGain(i)->Reset();
382 GetPedProfileHighGain(i)->Reset();
383 GetPedLEDRefProfileLowGain(i)->Reset();
384 GetPedLEDRefProfileHighGain(i)->Reset();
385 GetPeakProfileLowGain(i)->Reset();
386 GetPeakProfileHighGain(i)->Reset();
387 GetPeakHighGainHisto(i)->Reset();
388 GetDeadMap(i)->Reset();
390 if (!fPedestalLowGainDiff.IsEmpty()) {
391 //This means that the comparison profiles have been created.
393 GetPedProfileLowGainDiff(i)->Reset();
394 GetPedProfileHighGainDiff(i)->Reset();
395 GetPedLEDRefProfileLowGainDiff(i)->Reset();
396 GetPedLEDRefProfileHighGainDiff(i)->Reset();
397 GetPeakProfileLowGainDiff(i)->Reset();
398 GetPeakProfileHighGainDiff(i)->Reset();
400 GetPedProfileLowGainRatio(i)->Reset();
401 GetPedProfileHighGainRatio(i)->Reset();
402 GetPedLEDRefProfileLowGainRatio(i)->Reset();
403 GetPedLEDRefProfileHighGainRatio(i)->Reset();
404 GetPeakProfileLowGainRatio(i)->Reset();
405 GetPeakProfileHighGainRatio(i)->Reset();
412 fResurrectedTowers = 0;
414 //To think about: should fReference be deleted too?... let's not do it this time, at least...
417 // Parameter/cut handling
418 //_____________________________________________________________________
419 void AliCaloCalibPedestal::SetParametersFromFile(const char *parameterFile)
421 // Note: this method is a bit more complicated than it really has to be
422 // - allowing for multiple entries per line, arbitrary order of the
423 // different variables etc. But I wanted to try and do this in as
424 // correct a C++ way as I could (as an exercise).
426 static const string delimitor("::");
428 // open, check input file
429 ifstream in( parameterFile );
431 printf("in AliCaloCalibPedestal::SetParametersFromFile - Using default/run_time parameters.\n");
438 while ((in.rdstate() & ios::failbit) == 0 ) {
440 // Read into the raw char array and then construct a string
441 // to do the searching
442 in.getline(readline, 1024);
443 istringstream s(readline);
445 while ( ( s.rdstate() & ios::failbit ) == 0 ) {
450 // check stream status
451 if( s.rdstate() & ios::failbit ) break;
453 // skip rest of line if comments found
454 if( keyValue.substr( 0, 2 ) == "//" ) break;
456 // look for "::" in keyValue pair
457 size_t position = keyValue.find( delimitor );
458 if( position == string::npos ) {
459 printf("wrong format for key::value pair: %s\n", keyValue.c_str());
462 // split keyValue pair
463 string key( keyValue.substr( 0, position ) );
464 string value( keyValue.substr( position+delimitor.size(),
465 keyValue.size()-delimitor.size() ) );
467 // check value does not contain a new delimitor
468 if( value.find( delimitor ) != string::npos ) {
469 printf("wrong format for key::value pair: %s\n", keyValue.c_str());
472 // debug: check key value pair
473 // printf("AliCaloCalibPedestal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
475 // if the key matches with something we expect, we assign the new value
476 istringstream iss(value);
477 // the comparison strings defined at the beginning of this method
478 if ( (key == "fFirstPedestalSample") || (key == "fLastPedestalSample") || (key == "fDeadThreshold") || (key == "fWarningThreshold") || (key == "fWarningFraction") || (key == "fHotSigma") ) {
479 printf("AliCaloCalibPedestal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
481 if (key == "fFirstPedestalSample") {
482 iss >> fFirstPedestalSample;
484 else if (key == "fLastPedestalSample") {
485 iss >> fLastPedestalSample;
487 else if (key == "fDeadThreshold") {
488 iss >> fDeadThreshold;
490 else if (key == "fWarningThreshold") {
491 iss >> fWarningThreshold;
493 else if (key == "fWarningFraction") {
494 iss >> fWarningFraction;
496 else if (key == "fHotSigma") {
510 //_____________________________________________________________________
511 void AliCaloCalibPedestal::WriteParametersToFile(const char *parameterFile)
513 //Write parameters in file.
515 static const string delimitor("::");
516 ofstream out( parameterFile );
517 out << "// " << parameterFile << endl;
518 out << "fFirstPedestalSample" << "::" << fFirstPedestalSample << endl;
519 out << "fLastPedestalSample" << "::" << fLastPedestalSample << endl;
520 out << "fDeadThreshold" << "::" << fDeadThreshold << endl;
521 out << "fWarningThreshold" << "::" << fWarningThreshold << endl;
522 out << "fWarningFraction" << "::" << fWarningFraction << endl;
523 out << "fHotSigma" << "::" << fHotSigma << endl;
529 //_____________________________________________________________________
530 Bool_t AliCaloCalibPedestal::AddInfo(const AliCaloCalibPedestal *ped)
532 // just do this for the basic histograms/profiles that get filled in ProcessEvent
533 // may not have data for all modules, but let's just Add everything..
534 for (int i = 0; i < fModules; i++) {
535 GetPedProfileLowGain(i)->Add( ped->GetPedProfileLowGain(i) );
536 GetPedProfileHighGain(i)->Add( ped->GetPedProfileHighGain(i) );
537 GetPeakProfileLowGain(i)->Add( ped->GetPeakProfileLowGain(i) );
538 GetPeakProfileHighGain(i)->Add( ped->GetPeakProfileHighGain(i) );
539 GetPeakHighGainHisto(i)->Add( ped->GetPeakHighGainHisto(i) );
543 // DeadMap; Diff profiles etc would need to be redone after this operation
545 return kTRUE;//We succesfully added info from the supplied object
548 //_____________________________________________________________________
549 Bool_t AliCaloCalibPedestal::ProcessEvent(AliRawReader *rawReader)
551 // if fMapping is NULL the rawstream will crate its own mapping
552 AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
553 if (fDetType == kEmCal) {
554 rawReader->Select("EMCAL", 0, AliEMCALGeoParams::fgkLastAltroDDL) ; //select EMCAL DDL range
556 return ProcessEvent(&rawStream);
559 //_____________________________________________________________________
560 Bool_t AliCaloCalibPedestal::ProcessEvent(AliCaloRawStreamV3 *in)
562 // Method to process=analyze one event in the data stream
563 if (!in) return kFALSE; //Return right away if there's a null pointer
564 fNEvents++; // one more event
566 // indices for the reading
569 int i = 0; // sample counter
572 // start loop over input stream
573 while (in->NextDDL()) {
574 while (in->NextChannel()) {
577 int max = AliEMCALGeoParams::fgkSampleMin, min = AliEMCALGeoParams::fgkSampleMax; // min and max sample values
582 vector<int> pedSamples;
584 while (in->NextBunch()) {
585 const UShort_t *sig = in->GetSignals();
586 startBin = in->GetStartTimeBin();
587 nsamples += in->GetBunchLength();
588 for (i = 0; i < in->GetBunchLength(); i++) {
592 // check if it's a min or max value
593 if (sample < min) min = sample;
594 if (sample > max) max = sample;
596 // should we add it for the pedestal calculation?
597 if ( (fFirstPedestalSample<=time && time<=fLastPedestalSample) || // sample time in range
598 !fSelectPedestalSamples ) { // or we don't restrict the sample range.. - then we'll take all
599 pedSamples.push_back( sig[i] );
603 } // loop over samples in bunch
604 } // loop over bunches
606 if (nsamples > 0) { // this check is needed for when we have zero-supp. on, but not sparse readout
608 // it should be enough to check the SuperModule info for each DDL really, but let's keep it here for now
609 int arrayPos = in->GetModule(); //The modules are numbered starting from 0
610 if (arrayPos >= fModules) {
611 //TODO: return an error message, if appopriate (perhaps if debug>0?)
615 if (arrayPos < 0 || arrayPos >= fModules) {
616 printf("Oh no: arrayPos = %i.\n", arrayPos);
619 fNChanFills++; // one more channel found, and profile to be filled
620 //NOTE: coordinates are (column, row) for the profiles
621 if ( in->IsLowGain() ) {
622 //fill the low gain histograms
623 ((TProfile2D*)fPeakMinusPedLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min);
624 if (nPed>0) { // only fill pedestal info in case it could be calculated
625 for ( i=0; i<nPed; i++) {
626 ((TProfile2D*)fPedestalLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), pedSamples[i]);
630 else if ( in->IsHighGain() ) {
631 //fill the high gain ones
632 ((TProfile2D*)fPeakMinusPedHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min);
633 if (nPed>0) { // only fill pedestal info in case it could be calculated
634 for ( i=0; i<nPed; i++) {
635 ((TProfile2D*)fPedestalHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), pedSamples[i]);
638 // for warning checks
639 int idx = in->GetRow() + fRows * in->GetColumn();
640 ((TH2F*)fPeakMinusPedHighGainHisto[arrayPos])->Fill(idx, max - min);
642 else if ( in->IsLEDMonData() ) {
643 // for LED Mon data, the mapping class holds the gain info in the Row variable
644 // and the Strip number in the Column..
645 int gain = in->GetRow();
646 int stripId = in->GetColumn();
647 if (nPed>0 && stripId<fLEDRefs) {
649 for ( i=0; i<nPed; i++) {
650 ((TProfile*)fPedestalLEDRefLowGain[arrayPos])->Fill(stripId, pedSamples[i]);
654 for ( i=0; i<nPed; i++) {
655 ((TProfile*)fPedestalLEDRefHighGain[arrayPos])->Fill(stripId, pedSamples[i]);
661 } // nsamples>0 check, some data found for this channel; not only trailer/header
662 }// end while over channel
663 }//end while over DDL's, of input stream
665 in->Reset(); // just in case the next customer forgets to check if the stream was reset..
670 //_____________________________________________________________________
671 Bool_t AliCaloCalibPedestal::SaveHistograms(TString fileName, Bool_t saveEmptyHistos)
673 //Saves all the histograms (or profiles, to be accurate) to the designated file
675 TFile destFile(fileName, "recreate");
677 if (destFile.IsZombie()) {
683 for (int i = 0; i < fModules; i++) {
684 if( ((TProfile2D *)fPeakMinusPedLowGain[i])->GetEntries() || saveEmptyHistos) {
685 fPeakMinusPedLowGain[i]->Write();
687 if( ((TProfile2D *)fPeakMinusPedHighGain[i])->GetEntries() || saveEmptyHistos) {
688 fPeakMinusPedHighGain[i]->Write();
690 if( ((TProfile2D *)fPedestalLowGain[i])->GetEntries() || saveEmptyHistos) {
691 fPedestalLowGain[i]->Write();
693 if( ((TProfile2D *)fPedestalHighGain[i])->GetEntries() || saveEmptyHistos) {
694 fPedestalHighGain[i]->Write();
696 if( ((TProfile *)fPedestalLEDRefLowGain[i])->GetEntries() || saveEmptyHistos) {
697 fPedestalLEDRefLowGain[i]->Write();
699 if( ((TProfile *)fPedestalLEDRefHighGain[i])->GetEntries() || saveEmptyHistos) {
700 fPedestalLEDRefHighGain[i]->Write();
702 if( ((TH2F *)fPeakMinusPedHighGainHisto[i])->GetEntries() || saveEmptyHistos) {
703 fPeakMinusPedHighGainHisto[i]->Write();
713 //_____________________________________________________________________
714 Bool_t AliCaloCalibPedestal::LoadReferenceCalib(TString fileName, TString objectName)
717 //Make sure that the histograms created when loading the object are not destroyed as the file object is destroyed
718 TH1::AddDirectory(kFALSE);
720 TFile *sourceFile = new TFile(fileName);
721 if (sourceFile->IsZombie()) {
722 return kFALSE;//We couldn't load the reference
725 if (fReference) delete fReference;//Delete the reference object, if it already exists
728 fReference = (AliCaloCalibPedestal*)sourceFile->Get(objectName);
730 if (!fReference || !(fReference->InheritsFrom(AliCaloCalibPedestal::Class())) || (fReference->GetDetectorType() != fDetType)) {
731 if (fReference) delete fReference;//Delete the object, in case we had an object of the wrong type
738 //Reset the histogram ownership behaviour. NOTE: a better workaround would be good, since this may accidentally set AddDirectory to true, even
739 //if we are called by someone who has set it to false...
740 TH1::AddDirectory(kTRUE);
742 return kTRUE;//We succesfully loaded the object
746 //_____________________________________________________________________
747 Bool_t AliCaloCalibPedestal::SetReference(AliCaloCalibPedestal *ref)
749 if (fReference) delete fReference;//Delete the reference object, if it already exists
754 if (!fReference || (fReference->GetDetectorType() != fDetType)) {
755 if (fReference) delete fReference;//Delete the object, in case we had an object of the wrong type
760 return kTRUE;//We succesfully loaded the object
763 //_____________________________________________________________________
764 void AliCaloCalibPedestal::ValidateComparisonProfiles()
766 //Make sure the comparison histos exist
767 if (!fPedestalLowGainDiff.IsEmpty()) return; //The profiles already exist. We just check one, because they're all created at
771 //Then, loop for the requested number of modules
773 for (int i = 0; i < fModules; i++) {
774 //Pedestals, low gain
775 name = "hPedlowgainDiff";
777 title = "Pedestals difference, low gain, module ";
779 fPedestalLowGainDiff.Add(new TProfile2D(name, title,
780 fColumns, 0.0, fColumns,
781 fRows, fRowMin, fRowMax,"s"));
783 //Pedestals, high gain
784 name = "hPedhighgainDiff";
786 title = "Pedestals difference, high gain, module ";
788 fPedestalHighGainDiff.Add(new TProfile2D(name, title,
789 fColumns, 0.0, fColumns,
790 fRows, fRowMin, fRowMax,"s"));
792 //LED Ref/Mon pedestals, low gain
793 name = "hPedestalLEDReflowgainDiff";
795 title = "Pedestal difference LEDRef, low gain, module ";
797 fPedestalLEDRefLowGainDiff.Add(new TProfile(name, title,
798 fLEDRefs, 0.0, fLEDRefs, "s"));
800 //LED Ref/Mon pedestals, high gain
801 name = "hPedestalLEDRefhighgainDiff";
803 title = "Pedestal difference LEDRef, high gain, module ";
805 fPedestalLEDRefHighGainDiff.Add(new TProfile(name, title,
806 fLEDRefs, 0.0, fLEDRefs, "s"));
808 //Peak-Pedestals, high gain
809 name = "hPeakMinusPedhighgainDiff";
811 title = "Peak-Pedestal difference, high gain, module ";
813 fPeakMinusPedHighGainDiff.Add(new TProfile2D(name, title,
814 fColumns, 0.0, fColumns,
815 fRows, fRowMin, fRowMax,"s"));
817 //Peak-Pedestals, low gain
818 name = "hPeakMinusPedlowgainDiff";
820 title = "Peak-Pedestal difference, low gain, module ";
822 fPeakMinusPedLowGainDiff.Add(new TProfile2D(name, title,
823 fColumns, 0.0, fColumns,
824 fRows, fRowMin, fRowMax,"s"));
826 //Pedestals, low gain
827 name = "hPedlowgainRatio";
829 title = "Pedestals ratio, low gain, module ";
831 fPedestalLowGainRatio.Add(new TProfile2D(name, title,
832 fColumns, 0.0, fColumns,
833 fRows, fRowMin, fRowMax,"s"));
835 //Pedestals, high gain
836 name = "hPedhighgainRatio";
838 title = "Pedestals ratio, high gain, module ";
840 fPedestalHighGainRatio.Add(new TProfile2D(name, title,
841 fColumns, 0.0, fColumns,
842 fRows, fRowMin, fRowMax,"s"));
844 //LED Ref/Mon pedestals, low gain
845 name = "hPedestalLEDReflowgainRatio";
847 title = "Pedestal ratio LEDRef, low gain, module ";
849 fPedestalLEDRefLowGainRatio.Add(new TProfile(name, title,
850 fLEDRefs, 0.0, fLEDRefs, "s"));
852 //LED Ref/Mon pedestals, high gain
853 name = "hPedestalLEDRefhighgainRatio";
855 title = "Pedestal ratio LEDRef, high gain, module ";
857 fPedestalLEDRefHighGainRatio.Add(new TProfile(name, title,
858 fLEDRefs, 0.0, fLEDRefs, "s"));
860 //Peak-Pedestals, low gain
861 name = "hPeakMinusPedlowgainRatio";
863 title = "Peak-Pedestal ratio, low gain, module ";
865 fPeakMinusPedLowGainRatio.Add(new TProfile2D(name, title,
866 fColumns, 0.0, fColumns,
867 fRows, fRowMin, fRowMax,"s"));
869 //Peak-Pedestals, high gain
870 name = "hPeakMinusPedhighgainRatio";
872 title = "Peak-Pedestal ratio, high gain, module ";
874 fPeakMinusPedHighGainRatio.Add(new TProfile2D(name, title,
875 fColumns, 0.0, fColumns,
876 fRows, fRowMin, fRowMax,"s"));
878 }//end for nModules create the histograms
881 //_____________________________________________________________________
882 void AliCaloCalibPedestal::ComputeDiffAndRatio()
884 // calculate differences and ratios relative to a reference
885 ValidateComparisonProfiles();//Make sure the comparison histos exist
888 return;//Return if the reference object isn't loaded
894 for (int i = 0; i < fModules; i++) {
895 //For computing the difference, we cannot simply do TProfile2D->Add(), because that subtracts the sum of all entries,
896 //which means that the mean of the new profile will not be the difference of the means. So do it by hand:
897 for (int j = 0; j < fColumns; j++) {
898 for (int k = 0; k < fRows; k++) {
899 bin = ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->GetBin(j+1, k+1);//Note that we assume here that all histos have the same structure...
901 if (fReference->GetPeakProfileHighGain(i)->GetBinContent(bin) > 0) {
902 diff = GetPeakProfileHighGain(i)->GetBinContent(bin) - fReference->GetPeakProfileHighGain(i)->GetBinContent(bin);
903 ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinContent(bin, diff);
904 ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinEntries(bin, 1);
905 ratio = GetPeakProfileHighGain(i)->GetBinContent(bin) / fReference->GetPeakProfileHighGain(i)->GetBinContent(bin);
906 ((TProfile2D*)fPeakMinusPedHighGainRatio[i])->SetBinContent(bin, ratio);
907 ((TProfile2D*)fPeakMinusPedHighGainRatio[i])->SetBinEntries(bin, 1);
910 if (fReference->GetPeakProfileLowGain(i)->GetBinContent(bin) > 0) {
911 diff = GetPeakProfileLowGain(i)->GetBinContent(bin) - fReference->GetPeakProfileLowGain(i)->GetBinContent(bin);
912 ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinContent(bin, diff);
913 ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinEntries(bin, 1);
914 ratio = GetPeakProfileLowGain(i)->GetBinContent(bin) / fReference->GetPeakProfileLowGain(i)->GetBinContent(bin);
915 ((TProfile2D*)fPeakMinusPedLowGainRatio[i])->SetBinContent(bin, ratio);
916 ((TProfile2D*)fPeakMinusPedLowGainRatio[i])->SetBinEntries(bin, 1);
919 if (fReference->GetPedProfileHighGain(i)->GetBinContent(bin) > 0) {
920 diff = GetPedProfileHighGain(i)->GetBinContent(bin) - fReference->GetPedProfileHighGain(i)->GetBinContent(bin);
921 ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinContent(bin, diff);
922 ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinEntries(bin, 1);
923 ratio = GetPedProfileHighGain(i)->GetBinContent(bin) / fReference->GetPedProfileHighGain(i)->GetBinContent(bin);
924 ((TProfile2D*)fPedestalHighGainRatio[i])->SetBinContent(bin, ratio);
925 ((TProfile2D*)fPedestalHighGainRatio[i])->SetBinEntries(bin, 1);
928 if (fReference->GetPedProfileLowGain(i)->GetBinContent(bin) > 0) {
929 diff = GetPedProfileLowGain(i)->GetBinContent(bin) - fReference->GetPedProfileLowGain(i)->GetBinContent(bin);
930 ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinContent(bin, diff);
931 ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinEntries(bin, 1);
932 ratio = GetPedProfileLowGain(i)->GetBinContent(bin) / fReference->GetPedProfileLowGain(i)->GetBinContent(bin);
933 ((TProfile2D*)fPedestalLowGainRatio[i])->SetBinContent(bin, ratio);
934 ((TProfile2D*)fPedestalLowGainRatio[i])->SetBinEntries(bin, 1);
940 // same for LED Ref/Mon channels
941 for (int j = 0; j <= fLEDRefs; j++) {
942 bin = j+1;//Note that we assume here that all histos have the same structure...
944 if (fReference->GetPedLEDRefProfileHighGain(i)->GetBinContent(bin) > 0) {
945 diff = GetPedLEDRefProfileHighGain(i)->GetBinContent(bin) - fReference->GetPedLEDRefProfileHighGain(i)->GetBinContent(bin);
946 ((TProfile*)fPedestalLEDRefHighGainDiff[i])->SetBinContent(bin, diff);
947 ((TProfile*)fPedestalLEDRefHighGainDiff[i])->SetBinEntries(bin, 1);
948 ratio = GetPedLEDRefProfileHighGain(i)->GetBinContent(bin) / fReference->GetPedLEDRefProfileHighGain(i)->GetBinContent(bin);
949 ((TProfile*)fPedestalLEDRefHighGainRatio[i])->SetBinContent(bin, ratio);
950 ((TProfile*)fPedestalLEDRefHighGainRatio[i])->SetBinEntries(bin, 1);
953 if (fReference->GetPedLEDRefProfileLowGain(i)->GetBinContent(bin) > 0) {
954 diff = GetPedLEDRefProfileLowGain(i)->GetBinContent(bin) - fReference->GetPedLEDRefProfileLowGain(i)->GetBinContent(bin);
955 ((TProfile*)fPedestalLEDRefLowGainDiff[i])->SetBinContent(bin, diff);
956 ((TProfile*)fPedestalLEDRefLowGainDiff[i])->SetBinEntries(bin, 1);
957 ratio = GetPedLEDRefProfileLowGain(i)->GetBinContent(bin) / fReference->GetPedLEDRefProfileLowGain(i)->GetBinContent(bin);
958 ((TProfile*)fPedestalLEDRefLowGainRatio[i])->SetBinContent(bin, ratio);
959 ((TProfile*)fPedestalLEDRefLowGainRatio[i])->SetBinEntries(bin, 1);
968 //_____________________________________________________________________
969 void AliCaloCalibPedestal::ComputeHotAndWarningTowers(const char * hotMapFile)
970 { // look for hot/noisy towers
972 char name[512];//Quite a long temp buffer, just in case the filename includes a path
975 snprintf(name, 512, "%s.txt", hotMapFile);
976 fout = new ofstream(name);
977 if (!fout->is_open()) {
979 fout = 0;//Set the pointer to empty if the file was not opened
983 for(int i = 0; i < fModules; i++){
985 //first we compute the peak-pedestal distribution for each supermodule...
986 if( GetPeakHighGainHisto(i)->GetEntries() > 0 ) {
987 double min = GetPeakProfileHighGain(i)->GetBinContent(GetPeakProfileHighGain(i)->GetMinimumBin());
988 double max = GetPeakProfileHighGain(i)->GetBinContent(GetPeakProfileHighGain(i)->GetMaximumBin());
989 TH1D *hPeakFit = new TH1D(Form("hFit_%d", i), Form("hFit_%d", i), (int)((max-min)*10), min-1, max+1);
991 for (int j = 1; j <= fColumns; j++) {
992 for (int k = 1; k <= fRows; k++) {
993 hPeakFit->Fill(GetPeakProfileHighGain(i)->GetBinContent(j, k));
997 //...and then we try to fit it
998 double mean = hPeakFit->GetMean();
999 double sigma = hPeakFit->GetRMS();
1001 hPeakFit->Fit("gaus", "OQ", "", mean - 3*sigma, mean + 3*sigma);
1002 mean = hPeakFit->GetFunction("gaus")->GetParameter(1);
1003 sigma = hPeakFit->GetFunction("gaus")->GetParameter(2);
1005 catch (const std::exception & e) {
1006 printf("AliCaloCalibPedestal: TH1D PeakFit exception %s", e.what());
1012 //Then we look for warm/hot towers
1013 TH2F * hPeak2D = GetPeakHighGainHisto(i);
1014 hPeak2D->GetYaxis()->SetRangeUser( fWarningThreshold, hPeak2D->GetYaxis()->GetBinUpEdge(hPeak2D->GetNbinsY()) );
1017 int warnCounter = 0;
1018 for (int j = 1; j <= fColumns; j++) {
1019 for (int k = 1; k <= fRows; k++) {
1020 //we start looking for warm/warning towers...
1021 // histogram x-axis index
1022 idx = k-1 + fRows*(j-1); // this is what is used in the Fill call
1023 hPeak2D->GetXaxis()->SetRangeUser(idx, idx);
1024 warnCounter = (int) hPeak2D->Integral();
1025 if(warnCounter > fNEvents * fWarningFraction) {
1026 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kWarning);
1027 /* printf("mod %d col %d row %d warnCounter %d - status %d\n",
1028 i, j-1, k-1, warnCounter, (int) (kWarning)); */
1030 //...then we look for hot ones (towers whose values are greater than mean + X*sigma)
1031 if(GetPeakProfileHighGain(i)->GetBinContent(j, k) > mean + fHotSigma*sigma ) {
1032 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kHot);
1033 /* printf("mod %d col %d row %d binc %d - status %d\n",
1034 i, j-1, k-1, (int)(GetPeakProfileHighGain(i)->GetBinContent(j, k)), (int) (kHot)); */
1037 //Write the status to the hot/warm map file, if the file is open.
1038 // module - column - row - status (1=dead, 2= warm/warning , 3 = hot, see .h file enum)
1039 if (fout && ((TH2D*)fDeadMap[i])->GetBinContent(j, k) > 1) {
1044 << ((TH2D*)fDeadMap[i])->GetBinContent(j, k) << endl; }
1054 //_____________________________________________________________________
1055 void AliCaloCalibPedestal::ComputeDeadTowers(const char * deadMapFile)
1057 //Computes the number of dead towers etc etc into memory, after this you can call the GetDead... -functions
1061 ofstream * fout = 0;
1062 ofstream * diff = 0;
1063 char name[512];//Quite a long temp buffer, just in case the filename includes a path
1066 snprintf(name, 512, "%s.txt", deadMapFile);
1067 fout = new ofstream(name);
1068 snprintf(name, 512, "%sdiff.txt", deadMapFile);
1069 diff = new ofstream(name);
1070 if (!fout->is_open()) {
1072 fout = 0;//Set the pointer to empty if the file was not opened
1074 if (!diff->is_open()) {
1076 fout = 0;//Set the pointer to empty if the file was not opened
1080 for (int i = 0; i < fModules; i++) {
1081 if (GetPeakProfileHighGain(i)->GetEntries() > 0) { //don't care about empty histos
1082 for (int j = 1; j <= fColumns; j++) {
1083 for (int k = 1; k <= fRows; k++) {
1085 if (GetPeakProfileHighGain(i)->GetBinContent(j, k) < fDeadThreshold) {//It's dead
1086 countTot++;//One more dead total
1092 << "0" << endl;//Write the status to the deadmap file, if the file is open.
1095 if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) >= fDeadThreshold) {
1096 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kRecentlyDeceased);
1097 countNew++;//This tower wasn't dead before!
1099 ( *diff) << i << " "
1103 << "0" << endl;//Write the status to the deadmap difference file, if the file is open.
1107 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kDead);//This has been dead before. Nothing new
1110 else { //It's ALIVE!!
1111 //Don't bother with writing the live ones.
1112 if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) < fDeadThreshold) {
1113 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kResurrected);
1114 countRes++; //This tower was dead before => it's a miracle! :P
1120 << "1" << endl;//Write the status to the deadmap difference file, if the file is open.
1124 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kAlive);
1129 }//end for j/columns
1130 }//end if GetEntries >= 0
1139 fDeadTowers = countTot;
1140 fNewDeadTowers = countNew;
1141 fResurrectedTowers = countRes;
1144 //_____________________________________________________________________
1145 Bool_t AliCaloCalibPedestal::IsBadChannel(int imod, int icol, int irow) const
1147 //Check if channel is dead or hot.
1148 Int_t status = (Int_t) ( ((TH2D*)fDeadMap[imod])->GetBinContent(icol,irow) );
1149 if(status == kAlive)
1156 //_____________________________________________________________________
1157 void AliCaloCalibPedestal::SetChannelStatus(int imod, int icol, int irow, int status)
1159 //Set status of channel dead, hot, alive ...
1160 ((TH2D*)fDeadMap[imod])->SetBinContent(icol, irow, status);