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 // ValidateProfiles(); // not to be done in ctor; info from Axel N.
134 //_____________________________________________________________________
135 void AliCaloCalibPedestal::ValidateProfiles()
137 //Make sure the basic histos exist
138 if (!fPedestalLowGain.IsEmpty()) return; //The profiles already exist. We just check one, because they're all created at
141 //Then, loop for the requested number of modules
143 for (int i = 0; i < fModules; i++) {
144 //Pedestals, low gain
145 name = "hPedlowgain";
147 title = "Pedestals, low gain, module ";
149 fPedestalLowGain.Add(new TProfile2D(name, title,
150 fColumns, 0.0, fColumns,
151 fRows, fRowMin, fRowMax,"s"));
153 //Pedestals, high gain
154 name = "hPedhighgain";
156 title = "Pedestals, high gain, module ";
158 fPedestalHighGain.Add(new TProfile2D(name, title,
159 fColumns, 0.0, fColumns,
160 fRows, fRowMin, fRowMax,"s"));
162 //LED Ref/Mon pedestals, low gain
163 name = "hPedestalLEDReflowgain";
165 title = "Pedestal LEDRef, low gain, module ";
167 fPedestalLEDRefLowGain.Add(new TProfile(name, title,
168 fLEDRefs, 0.0, fLEDRefs, "s"));
170 //LED Ref/Mon pedestals, high gain
171 name = "hPedestalLEDRefhighgain";
173 title = "Pedestal LEDRef, high gain, module ";
175 fPedestalLEDRefHighGain.Add(new TProfile(name, title,
176 fLEDRefs, 0.0, fLEDRefs, "s"));
178 //Peak-Pedestals, low gain
179 name = "hPeakMinusPedlowgain";
181 title = "Peak-Pedestal, low gain, module ";
183 fPeakMinusPedLowGain.Add(new TProfile2D(name, title,
184 fColumns, 0.0, fColumns,
185 fRows, fRowMin, fRowMax,"s"));
187 //Peak-Pedestals, high gain
188 name = "hPeakMinusPedhighgain";
190 title = "Peak-Pedestal, high gain, module ";
192 fPeakMinusPedHighGain.Add(new TProfile2D(name, title,
193 fColumns, 0.0, fColumns,
194 fRows, fRowMin, fRowMax,"s"));
196 //Peak-Pedestals, high gain - TH2F histo
197 name = "hPeakMinusPedhighgainHisto";
199 title = "Peak-Pedestal, high gain, module ";
201 fPeakMinusPedHighGainHisto.Add(new TH2F(name, title,
202 fColumns*fRows, 0.0, fColumns*fRows,
207 title = "Dead map, module ";
209 fDeadMap.Add(new TH2D(name, title, fColumns, 0.0, fColumns,
210 fRows, fRowMin, fRowMax));
212 }//end for nModules create the histograms
214 CompressAndSetOwner();
217 //_____________________________________________________________________
218 void AliCaloCalibPedestal::CompressAndSetOwner()
220 //Compress the arrays, in order to remove the empty objects (a 16 slot array is created by default)
221 fPedestalLowGain.Compress();
222 fPedestalHighGain.Compress();
223 fPedestalLEDRefLowGain.Compress();
224 fPedestalLEDRefHighGain.Compress();
225 fPeakMinusPedLowGain.Compress();
226 fPeakMinusPedHighGain.Compress();
227 fPeakMinusPedHighGainHisto.Compress();
230 // set owner ship for everyone
231 fPedestalLowGain.SetOwner(kTRUE);
232 fPedestalHighGain.SetOwner(kTRUE);
233 fPedestalLEDRefLowGain.SetOwner(kTRUE);
234 fPedestalLEDRefHighGain.SetOwner(kTRUE);
235 fPeakMinusPedLowGain.SetOwner(kTRUE);
236 fPeakMinusPedHighGain.SetOwner(kTRUE);
237 fPeakMinusPedHighGainHisto.SetOwner(kTRUE);
238 fPedestalLowGainDiff.SetOwner(kTRUE);
239 fPedestalHighGainDiff.SetOwner(kTRUE);
240 fPedestalLEDRefLowGainDiff.SetOwner(kTRUE);
241 fPedestalLEDRefHighGainDiff.SetOwner(kTRUE);
242 fPeakMinusPedLowGainDiff.SetOwner(kTRUE);
243 fPeakMinusPedHighGainDiff.SetOwner(kTRUE);
244 fPedestalLowGainRatio.SetOwner(kTRUE);
245 fPedestalHighGainRatio.SetOwner(kTRUE);
246 fPedestalLEDRefLowGainRatio.SetOwner(kTRUE);
247 fPedestalLEDRefHighGainRatio.SetOwner(kTRUE);
248 fPeakMinusPedLowGainRatio.SetOwner(kTRUE);
249 fPeakMinusPedHighGainRatio.SetOwner(kTRUE);
250 fDeadMap.SetOwner(kTRUE);
254 //_____________________________________________________________________
255 AliCaloCalibPedestal::~AliCaloCalibPedestal()
259 if (fReference) delete fReference;//Delete the reference object, if it has been loaded
261 // delete also TObjArray's
262 fPedestalLowGain.Delete();
263 fPedestalHighGain.Delete();
264 fPedestalLEDRefLowGain.Delete();
265 fPedestalLEDRefHighGain.Delete();
266 fPeakMinusPedLowGain.Delete();
267 fPeakMinusPedHighGain.Delete();
268 fPeakMinusPedHighGainHisto.Delete();
269 fPedestalLowGainDiff.Delete();
270 fPedestalHighGainDiff.Delete();
271 fPedestalLEDRefLowGainDiff.Delete();
272 fPedestalLEDRefHighGainDiff.Delete();
273 fPeakMinusPedLowGainDiff.Delete();
274 fPeakMinusPedHighGainDiff.Delete();
275 fPedestalLowGainRatio.Delete();
276 fPedestalHighGainRatio.Delete();
277 fPedestalLEDRefLowGainRatio.Delete();
278 fPedestalLEDRefHighGainRatio.Delete();
279 fPeakMinusPedLowGainRatio.Delete();
280 fPeakMinusPedHighGainRatio.Delete();
286 //_____________________________________________________________________
287 AliCaloCalibPedestal::AliCaloCalibPedestal(AliCaloCalibPedestal &ped) :
291 fPedestalLEDRefLowGain(),
292 fPedestalLEDRefHighGain(),
293 fPeakMinusPedLowGain(),
294 fPeakMinusPedHighGain(),
295 fPeakMinusPedHighGainHisto(),
296 fPedestalLowGainDiff(),
297 fPedestalHighGainDiff(),
298 fPedestalLEDRefLowGainDiff(),
299 fPedestalLEDRefHighGainDiff(),
300 fPeakMinusPedLowGainDiff(),
301 fPeakMinusPedHighGainDiff(),
302 fPedestalLowGainRatio(),
303 fPedestalHighGainRatio(),
304 fPedestalLEDRefLowGainRatio(),
305 fPedestalLEDRefHighGainRatio(),
306 fPeakMinusPedLowGainRatio(),
307 fPeakMinusPedHighGainRatio(),
309 fNEvents(ped.GetNEvents()),
310 fNChanFills(ped.GetNChanFills()),
311 fDeadTowers(ped.GetDeadTowerCount()),
312 fNewDeadTowers(ped.GetDeadTowerNew()),
313 fResurrectedTowers(ped.GetDeadTowerResurrected()),
314 fReference( 0 ), //! note that we do not try to copy the reference info here
315 fDetType(ped.GetDetectorType()),
316 fColumns(ped.GetColumns()),
317 fRows(ped.GetRows()),
318 fLEDRefs(ped.GetLEDRefs()),
319 fModules(ped.GetModules()),
320 fRowMin(ped.GetRowMin()),
321 fRowMax(ped.GetRowMax()),
322 fRowMultiplier(ped.GetRowMultiplier()),
323 fCaloString(ped.GetCaloString()),
324 fMapping(NULL), //! note that we are not copying the map info
325 fRunNumber(ped.GetRunNumber()),
326 fSelectPedestalSamples(ped.GetSelectPedestalSamples()),
327 fFirstPedestalSample(ped.GetFirstPedestalSample()),
328 fLastPedestalSample(ped.GetLastPedestalSample()),
329 fDeadThreshold(ped.GetDeadThreshold()),
330 fWarningThreshold(ped.GetWarningThreshold()),
331 fWarningFraction(ped.GetWarningFraction()),
332 fHotSigma(ped.GetHotSigma())
334 // Then the ObjArray ones; we add the histograms rather than trying TObjArray = assignment
335 //DS: this has not really been tested yet..
336 for (int i = 0; i < fModules; i++) {
337 fPedestalLowGain.Add( ped.GetPedProfileLowGain(i) );
338 fPedestalHighGain.Add( ped.GetPedProfileHighGain(i) );
339 fPedestalLEDRefLowGain.Add( ped.GetPedLEDRefProfileLowGain(i) );
340 fPedestalLEDRefHighGain.Add( ped.GetPedLEDRefProfileHighGain(i) );
341 fPeakMinusPedLowGain.Add( ped.GetPeakProfileLowGain(i) );
342 fPeakMinusPedHighGain.Add( ped.GetPeakProfileHighGain(i) );
343 fPeakMinusPedHighGainHisto.Add( ped.GetPeakHighGainHisto(i) );
345 fDeadMap.Add( ped.GetDeadMap(i) );
348 CompressAndSetOwner();
351 // assignment operator; use copy ctor to make life easy..
352 //_____________________________________________________________________
353 AliCaloCalibPedestal& AliCaloCalibPedestal::operator = (AliCaloCalibPedestal &source)
355 // assignment operator; use copy ctor
356 if (&source == this) return *this;
358 new (this) AliCaloCalibPedestal(source);
362 //_____________________________________________________________________
363 void AliCaloCalibPedestal::Reset()
364 { // Reset all arrays/histograms
365 ValidateProfiles(); // make sure histos/profiles exist
366 for (int i = 0; i < fModules; i++) {
367 GetPedProfileLowGain(i)->Reset();
368 GetPedProfileHighGain(i)->Reset();
369 GetPedLEDRefProfileLowGain(i)->Reset();
370 GetPedLEDRefProfileHighGain(i)->Reset();
371 GetPeakProfileLowGain(i)->Reset();
372 GetPeakProfileHighGain(i)->Reset();
373 GetPeakHighGainHisto(i)->Reset();
374 GetDeadMap(i)->Reset();
376 if (!fPedestalLowGainDiff.IsEmpty()) {
377 //This means that the comparison profiles have been created.
379 GetPedProfileLowGainDiff(i)->Reset();
380 GetPedProfileHighGainDiff(i)->Reset();
381 GetPedLEDRefProfileLowGainDiff(i)->Reset();
382 GetPedLEDRefProfileHighGainDiff(i)->Reset();
383 GetPeakProfileLowGainDiff(i)->Reset();
384 GetPeakProfileHighGainDiff(i)->Reset();
386 GetPedProfileLowGainRatio(i)->Reset();
387 GetPedProfileHighGainRatio(i)->Reset();
388 GetPedLEDRefProfileLowGainRatio(i)->Reset();
389 GetPedLEDRefProfileHighGainRatio(i)->Reset();
390 GetPeakProfileLowGainRatio(i)->Reset();
391 GetPeakProfileHighGainRatio(i)->Reset();
398 fResurrectedTowers = 0;
400 //To think about: should fReference be deleted too?... let's not do it this time, at least...
403 // Parameter/cut handling
404 //_____________________________________________________________________
405 void AliCaloCalibPedestal::SetParametersFromFile(const char *parameterFile)
407 // Note: this method is a bit more complicated than it really has to be
408 // - allowing for multiple entries per line, arbitrary order of the
409 // different variables etc. But I wanted to try and do this in as
410 // correct a C++ way as I could (as an exercise).
412 static const string delimitor("::");
414 // open, check input file
415 ifstream in( parameterFile );
417 printf("in AliCaloCalibPedestal::SetParametersFromFile - Using default/run_time parameters.\n");
424 while ((in.rdstate() & ios::failbit) == 0 ) {
426 // Read into the raw char array and then construct a string
427 // to do the searching
428 in.getline(readline, 1024);
429 istringstream s(readline);
431 while ( ( s.rdstate() & ios::failbit ) == 0 ) {
436 // check stream status
437 if( ( s.rdstate() & ios::failbit ) == ios::failbit) break;
439 // skip rest of line if comments found
440 if( keyValue.substr( 0, 2 ) == "//" ) break;
442 // look for "::" in keyValue pair
443 size_t position = keyValue.find( delimitor );
444 if( position == string::npos ) {
445 printf("wrong format for key::value pair: %s\n", keyValue.c_str());
448 // split keyValue pair
449 string key( keyValue.substr( 0, position ) );
450 string value( keyValue.substr( position+delimitor.size(),
451 keyValue.size()-delimitor.size() ) );
453 // check value does not contain a new delimitor
454 if( value.find( delimitor ) != string::npos ) {
455 printf("wrong format for key::value pair: %s\n", keyValue.c_str());
458 // debug: check key value pair
459 // printf("AliCaloCalibPedestal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
461 // if the key matches with something we expect, we assign the new value
462 istringstream iss(value);
463 // the comparison strings defined at the beginning of this method
464 if ( (key == "fFirstPedestalSample") || (key == "fLastPedestalSample") || (key == "fDeadThreshold") || (key == "fWarningThreshold") || (key == "fWarningFraction") || (key == "fHotSigma") ) {
465 printf("AliCaloCalibPedestal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
467 if (key == "fFirstPedestalSample") {
468 iss >> fFirstPedestalSample;
470 else if (key == "fLastPedestalSample") {
471 iss >> fLastPedestalSample;
473 else if (key == "fDeadThreshold") {
474 iss >> fDeadThreshold;
476 else if (key == "fWarningThreshold") {
477 iss >> fWarningThreshold;
479 else if (key == "fWarningFraction") {
480 iss >> fWarningFraction;
482 else if (key == "fHotSigma") {
496 //_____________________________________________________________________
497 void AliCaloCalibPedestal::WriteParametersToFile(const char *parameterFile)
499 //Write parameters in file.
501 static const string delimitor("::");
502 ofstream out( parameterFile );
503 out << "// " << parameterFile << endl;
504 out << "fFirstPedestalSample" << "::" << fFirstPedestalSample << endl;
505 out << "fLastPedestalSample" << "::" << fLastPedestalSample << endl;
506 out << "fDeadThreshold" << "::" << fDeadThreshold << endl;
507 out << "fWarningThreshold" << "::" << fWarningThreshold << endl;
508 out << "fWarningFraction" << "::" << fWarningFraction << endl;
509 out << "fHotSigma" << "::" << fHotSigma << endl;
515 //_____________________________________________________________________
516 Bool_t AliCaloCalibPedestal::AddInfo(AliCaloCalibPedestal *ped)
518 // just do this for the basic histograms/profiles that get filled in ProcessEvent
519 // may not have data for all modules, but let's just Add everything..
520 ValidateProfiles(); // make sure histos/profiles exist
522 for (int i = 0; i < fModules; i++) {
523 GetPedProfileLowGain(i)->Add( ped->GetPedProfileLowGain(i) );
524 GetPedProfileHighGain(i)->Add( ped->GetPedProfileHighGain(i) );
525 GetPedLEDRefProfileLowGain(i)->Add( ped->GetPedLEDRefProfileLowGain(i) );
526 GetPedLEDRefProfileHighGain(i)->Add( ped->GetPedLEDRefProfileHighGain(i) );
527 GetPeakProfileLowGain(i)->Add( ped->GetPeakProfileLowGain(i) );
528 GetPeakProfileHighGain(i)->Add( ped->GetPeakProfileHighGain(i) );
529 GetPeakHighGainHisto(i)->Add( ped->GetPeakHighGainHisto(i) );
532 // We should also copy other pieces of info: counters and parameters
533 // (not number of columns and rows etc which should be the same)
534 // note that I just assign them here rather than Add them, but we
535 // normally just Add (e.g. in Preprocessor) one object so this should be fine.
536 fNEvents = ped->GetNEvents();
537 fNChanFills = ped->GetNChanFills();
538 fDeadTowers = ped->GetDeadTowerCount();
539 fNewDeadTowers = ped->GetDeadTowerNew();
540 fResurrectedTowers = ped->GetDeadTowerResurrected();
541 fRunNumber = ped->GetRunNumber();
542 fSelectPedestalSamples = ped->GetSelectPedestalSamples();
543 fFirstPedestalSample = ped->GetFirstPedestalSample();
544 fLastPedestalSample = ped->GetLastPedestalSample();
545 fDeadThreshold = ped->GetDeadThreshold();
546 fWarningThreshold = ped->GetWarningThreshold();
547 fWarningFraction = ped->GetWarningFraction();
548 fHotSigma = ped->GetHotSigma();
550 // DeadMap; Diff profiles etc would need to be redone after this operation
552 return kTRUE;//We succesfully added info from the supplied object
555 //_____________________________________________________________________
556 Bool_t AliCaloCalibPedestal::ProcessEvent(AliRawReader *rawReader)
558 // if fMapping is NULL the rawstream will crate its own mapping
559 AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
560 if (fDetType == kEmCal) {
561 rawReader->Select("EMCAL", 0, AliEMCALGeoParams::fgkLastAltroDDL) ; //select EMCAL DDL range
563 return ProcessEvent(&rawStream);
566 //_____________________________________________________________________
567 Bool_t AliCaloCalibPedestal::ProcessEvent(AliCaloRawStreamV3 *in)
569 // Method to process=analyze one event in the data stream
570 if (!in) return kFALSE; //Return right away if there's a null pointer
571 in->Reset(); // just in case the next customer forgets to check if the stream was reset..
573 fNEvents++; // one more event
575 if (fNEvents==1) ValidateProfiles(); // 1st event, make sure histos/profiles exist
577 // indices for the reading
580 int i = 0; // sample counter
583 // start loop over input stream
584 while (in->NextDDL()) {
585 while (in->NextChannel()) {
588 int max = AliEMCALGeoParams::fgkSampleMin, min = AliEMCALGeoParams::fgkSampleMax; // min and max sample values
593 vector<int> pedSamples;
595 while (in->NextBunch()) {
596 const UShort_t *sig = in->GetSignals();
597 startBin = in->GetStartTimeBin();
598 nsamples += in->GetBunchLength();
599 for (i = 0; i < in->GetBunchLength(); i++) {
603 // check if it's a min or max value
604 if (sample < min) min = sample;
605 if (sample > max) max = sample;
607 // should we add it for the pedestal calculation?
608 if ( (fFirstPedestalSample<=time && time<=fLastPedestalSample) || // sample time in range
609 !fSelectPedestalSamples ) { // or we don't restrict the sample range.. - then we'll take all
610 pedSamples.push_back( sig[i] );
614 } // loop over samples in bunch
615 } // loop over bunches
617 if (nsamples > 0) { // this check is needed for when we have zero-supp. on, but not sparse readout
619 // it should be enough to check the SuperModule info for each DDL really, but let's keep it here for now
620 int arrayPos = in->GetModule(); //The modules are numbered starting from 0
621 if (arrayPos >= fModules) {
622 //TODO: return an error message, if appopriate (perhaps if debug>0?)
626 if (arrayPos < 0 || arrayPos >= fModules) {
627 printf("Oh no: arrayPos = %i.\n", arrayPos);
630 fNChanFills++; // one more channel found, and profile to be filled
631 //NOTE: coordinates are (column, row) for the profiles
632 if ( in->IsLowGain() ) {
633 //fill the low gain histograms
634 ((TProfile2D*)fPeakMinusPedLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min);
635 if (nPed>0) { // only fill pedestal info in case it could be calculated
636 for ( i=0; i<nPed; i++) {
637 ((TProfile2D*)fPedestalLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), pedSamples[i]);
641 else if ( in->IsHighGain() ) {
642 //fill the high gain ones
643 ((TProfile2D*)fPeakMinusPedHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min);
644 if (nPed>0) { // only fill pedestal info in case it could be calculated
645 for ( i=0; i<nPed; i++) {
646 ((TProfile2D*)fPedestalHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), pedSamples[i]);
649 // for warning checks
650 int idx = in->GetRow() + fRows * in->GetColumn();
651 ((TH2F*)fPeakMinusPedHighGainHisto[arrayPos])->Fill(idx, max - min);
653 else if ( in->IsLEDMonData() ) {
654 // for LED Mon data, the mapping class holds the gain info in the Row variable
655 // and the Strip number in the Column..
656 int gain = in->GetRow();
657 int stripId = in->GetColumn();
658 if (nPed>0 && stripId<fLEDRefs) {
660 for ( i=0; i<nPed; i++) {
661 ((TProfile*)fPedestalLEDRefLowGain[arrayPos])->Fill(stripId, pedSamples[i]);
665 for ( i=0; i<nPed; i++) {
666 ((TProfile*)fPedestalLEDRefHighGain[arrayPos])->Fill(stripId, pedSamples[i]);
672 } // nsamples>0 check, some data found for this channel; not only trailer/header
673 }// end while over channel
674 }//end while over DDL's, of input stream
680 //_____________________________________________________________________
681 Bool_t AliCaloCalibPedestal::SaveHistograms(TString fileName, Bool_t saveEmptyHistos)
683 //Saves all the histograms (or profiles, to be accurate) to the designated file
684 ValidateProfiles(); // make sure histos/profiles exist
685 TFile destFile(fileName, "recreate");
687 if (destFile.IsZombie()) {
693 for (int i = 0; i < fModules; i++) {
694 if( ((TProfile2D *)fPeakMinusPedLowGain[i])->GetEntries() || saveEmptyHistos) {
695 fPeakMinusPedLowGain[i]->Write();
697 if( ((TProfile2D *)fPeakMinusPedHighGain[i])->GetEntries() || saveEmptyHistos) {
698 fPeakMinusPedHighGain[i]->Write();
700 if( ((TProfile2D *)fPedestalLowGain[i])->GetEntries() || saveEmptyHistos) {
701 fPedestalLowGain[i]->Write();
703 if( ((TProfile2D *)fPedestalHighGain[i])->GetEntries() || saveEmptyHistos) {
704 fPedestalHighGain[i]->Write();
706 if( ((TProfile *)fPedestalLEDRefLowGain[i])->GetEntries() || saveEmptyHistos) {
707 fPedestalLEDRefLowGain[i]->Write();
709 if( ((TProfile *)fPedestalLEDRefHighGain[i])->GetEntries() || saveEmptyHistos) {
710 fPedestalLEDRefHighGain[i]->Write();
712 if( ((TH2F *)fPeakMinusPedHighGainHisto[i])->GetEntries() || saveEmptyHistos) {
713 fPeakMinusPedHighGainHisto[i]->Write();
723 //_____________________________________________________________________
724 Bool_t AliCaloCalibPedestal::LoadReferenceCalib(TString fileName, TString objectName)
727 //Make sure that the histograms created when loading the object are not destroyed as the file object is destroyed
728 TH1::AddDirectory(kFALSE);
730 TFile *sourceFile = new TFile(fileName);
731 if (sourceFile->IsZombie()) {
732 return kFALSE;//We couldn't load the reference
735 if (fReference) delete fReference;//Delete the reference object, if it already exists
738 fReference = (AliCaloCalibPedestal*)sourceFile->Get(objectName);
740 if (!fReference || !(fReference->InheritsFrom(AliCaloCalibPedestal::Class())) || (fReference->GetDetectorType() != fDetType)) {
741 if (fReference) delete fReference;//Delete the object, in case we had an object of the wrong type
748 //Reset the histogram ownership behaviour. NOTE: a better workaround would be good, since this may accidentally set AddDirectory to true, even
749 //if we are called by someone who has set it to false...
750 TH1::AddDirectory(kTRUE);
752 return kTRUE;//We succesfully loaded the object
756 //_____________________________________________________________________
757 Bool_t AliCaloCalibPedestal::SetReference(AliCaloCalibPedestal *ref)
758 { // set reference object
759 if (fReference) delete fReference;//Delete the reference object, if it already exists
764 if (!fReference || (fReference->GetDetectorType() != fDetType)) {
765 if (fReference) delete fReference;//Delete the object, in case we had an object of the wrong type
770 return kTRUE;//We succesfully loaded the object
773 //_____________________________________________________________________
774 void AliCaloCalibPedestal::ValidateComparisonProfiles()
776 //Make sure the comparison histos exist
777 if (!fPedestalLowGainDiff.IsEmpty()) return; //The profiles already exist. We just check one, because they're all created at
781 //Then, loop for the requested number of modules
783 for (int i = 0; i < fModules; i++) {
784 //Pedestals, low gain
785 name = "hPedlowgainDiff";
787 title = "Pedestals difference, low gain, module ";
789 fPedestalLowGainDiff.Add(new TProfile2D(name, title,
790 fColumns, 0.0, fColumns,
791 fRows, fRowMin, fRowMax,"s"));
793 //Pedestals, high gain
794 name = "hPedhighgainDiff";
796 title = "Pedestals difference, high gain, module ";
798 fPedestalHighGainDiff.Add(new TProfile2D(name, title,
799 fColumns, 0.0, fColumns,
800 fRows, fRowMin, fRowMax,"s"));
802 //LED Ref/Mon pedestals, low gain
803 name = "hPedestalLEDReflowgainDiff";
805 title = "Pedestal difference LEDRef, low gain, module ";
807 fPedestalLEDRefLowGainDiff.Add(new TProfile(name, title,
808 fLEDRefs, 0.0, fLEDRefs, "s"));
810 //LED Ref/Mon pedestals, high gain
811 name = "hPedestalLEDRefhighgainDiff";
813 title = "Pedestal difference LEDRef, high gain, module ";
815 fPedestalLEDRefHighGainDiff.Add(new TProfile(name, title,
816 fLEDRefs, 0.0, fLEDRefs, "s"));
818 //Peak-Pedestals, high gain
819 name = "hPeakMinusPedhighgainDiff";
821 title = "Peak-Pedestal difference, high gain, module ";
823 fPeakMinusPedHighGainDiff.Add(new TProfile2D(name, title,
824 fColumns, 0.0, fColumns,
825 fRows, fRowMin, fRowMax,"s"));
827 //Peak-Pedestals, low gain
828 name = "hPeakMinusPedlowgainDiff";
830 title = "Peak-Pedestal difference, low gain, module ";
832 fPeakMinusPedLowGainDiff.Add(new TProfile2D(name, title,
833 fColumns, 0.0, fColumns,
834 fRows, fRowMin, fRowMax,"s"));
836 //Pedestals, low gain
837 name = "hPedlowgainRatio";
839 title = "Pedestals ratio, low gain, module ";
841 fPedestalLowGainRatio.Add(new TProfile2D(name, title,
842 fColumns, 0.0, fColumns,
843 fRows, fRowMin, fRowMax,"s"));
845 //Pedestals, high gain
846 name = "hPedhighgainRatio";
848 title = "Pedestals ratio, high gain, module ";
850 fPedestalHighGainRatio.Add(new TProfile2D(name, title,
851 fColumns, 0.0, fColumns,
852 fRows, fRowMin, fRowMax,"s"));
854 //LED Ref/Mon pedestals, low gain
855 name = "hPedestalLEDReflowgainRatio";
857 title = "Pedestal ratio LEDRef, low gain, module ";
859 fPedestalLEDRefLowGainRatio.Add(new TProfile(name, title,
860 fLEDRefs, 0.0, fLEDRefs, "s"));
862 //LED Ref/Mon pedestals, high gain
863 name = "hPedestalLEDRefhighgainRatio";
865 title = "Pedestal ratio LEDRef, high gain, module ";
867 fPedestalLEDRefHighGainRatio.Add(new TProfile(name, title,
868 fLEDRefs, 0.0, fLEDRefs, "s"));
870 //Peak-Pedestals, low gain
871 name = "hPeakMinusPedlowgainRatio";
873 title = "Peak-Pedestal ratio, low gain, module ";
875 fPeakMinusPedLowGainRatio.Add(new TProfile2D(name, title,
876 fColumns, 0.0, fColumns,
877 fRows, fRowMin, fRowMax,"s"));
879 //Peak-Pedestals, high gain
880 name = "hPeakMinusPedhighgainRatio";
882 title = "Peak-Pedestal ratio, high gain, module ";
884 fPeakMinusPedHighGainRatio.Add(new TProfile2D(name, title,
885 fColumns, 0.0, fColumns,
886 fRows, fRowMin, fRowMax,"s"));
888 }//end for nModules create the histograms
891 //_____________________________________________________________________
892 void AliCaloCalibPedestal::ComputeDiffAndRatio()
893 { // calculate differences and ratios relative to a reference
894 ValidateProfiles(); // make sure histos/profiles exist
895 ValidateComparisonProfiles();//Make sure the comparison histos exist
898 return;//Return if the reference object isn't loaded
904 for (int i = 0; i < fModules; i++) {
905 //For computing the difference, we cannot simply do TProfile2D->Add(), because that subtracts the sum of all entries,
906 //which means that the mean of the new profile will not be the difference of the means. So do it by hand:
907 for (int j = 0; j < fColumns; j++) {
908 for (int k = 0; k < fRows; k++) {
909 bin = ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->GetBin(j+1, k+1);//Note that we assume here that all histos have the same structure...
911 if (fReference->GetPeakProfileHighGain(i)->GetBinContent(bin) > 0) {
912 diff = GetPeakProfileHighGain(i)->GetBinContent(bin) - fReference->GetPeakProfileHighGain(i)->GetBinContent(bin);
913 ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinContent(bin, diff);
914 ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinEntries(bin, 1);
915 ratio = GetPeakProfileHighGain(i)->GetBinContent(bin) / fReference->GetPeakProfileHighGain(i)->GetBinContent(bin);
916 ((TProfile2D*)fPeakMinusPedHighGainRatio[i])->SetBinContent(bin, ratio);
917 ((TProfile2D*)fPeakMinusPedHighGainRatio[i])->SetBinEntries(bin, 1);
920 if (fReference->GetPeakProfileLowGain(i)->GetBinContent(bin) > 0) {
921 diff = GetPeakProfileLowGain(i)->GetBinContent(bin) - fReference->GetPeakProfileLowGain(i)->GetBinContent(bin);
922 ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinContent(bin, diff);
923 ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinEntries(bin, 1);
924 ratio = GetPeakProfileLowGain(i)->GetBinContent(bin) / fReference->GetPeakProfileLowGain(i)->GetBinContent(bin);
925 ((TProfile2D*)fPeakMinusPedLowGainRatio[i])->SetBinContent(bin, ratio);
926 ((TProfile2D*)fPeakMinusPedLowGainRatio[i])->SetBinEntries(bin, 1);
929 if (fReference->GetPedProfileHighGain(i)->GetBinContent(bin) > 0) {
930 diff = GetPedProfileHighGain(i)->GetBinContent(bin) - fReference->GetPedProfileHighGain(i)->GetBinContent(bin);
931 ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinContent(bin, diff);
932 ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinEntries(bin, 1);
933 ratio = GetPedProfileHighGain(i)->GetBinContent(bin) / fReference->GetPedProfileHighGain(i)->GetBinContent(bin);
934 ((TProfile2D*)fPedestalHighGainRatio[i])->SetBinContent(bin, ratio);
935 ((TProfile2D*)fPedestalHighGainRatio[i])->SetBinEntries(bin, 1);
938 if (fReference->GetPedProfileLowGain(i)->GetBinContent(bin) > 0) {
939 diff = GetPedProfileLowGain(i)->GetBinContent(bin) - fReference->GetPedProfileLowGain(i)->GetBinContent(bin);
940 ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinContent(bin, diff);
941 ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinEntries(bin, 1);
942 ratio = GetPedProfileLowGain(i)->GetBinContent(bin) / fReference->GetPedProfileLowGain(i)->GetBinContent(bin);
943 ((TProfile2D*)fPedestalLowGainRatio[i])->SetBinContent(bin, ratio);
944 ((TProfile2D*)fPedestalLowGainRatio[i])->SetBinEntries(bin, 1);
950 // same for LED Ref/Mon channels
951 for (int j = 0; j <= fLEDRefs; j++) {
952 bin = j+1;//Note that we assume here that all histos have the same structure...
954 if (fReference->GetPedLEDRefProfileHighGain(i)->GetBinContent(bin) > 0) {
955 diff = GetPedLEDRefProfileHighGain(i)->GetBinContent(bin) - fReference->GetPedLEDRefProfileHighGain(i)->GetBinContent(bin);
956 ((TProfile*)fPedestalLEDRefHighGainDiff[i])->SetBinContent(bin, diff);
957 ((TProfile*)fPedestalLEDRefHighGainDiff[i])->SetBinEntries(bin, 1);
958 ratio = GetPedLEDRefProfileHighGain(i)->GetBinContent(bin) / fReference->GetPedLEDRefProfileHighGain(i)->GetBinContent(bin);
959 ((TProfile*)fPedestalLEDRefHighGainRatio[i])->SetBinContent(bin, ratio);
960 ((TProfile*)fPedestalLEDRefHighGainRatio[i])->SetBinEntries(bin, 1);
963 if (fReference->GetPedLEDRefProfileLowGain(i)->GetBinContent(bin) > 0) {
964 diff = GetPedLEDRefProfileLowGain(i)->GetBinContent(bin) - fReference->GetPedLEDRefProfileLowGain(i)->GetBinContent(bin);
965 ((TProfile*)fPedestalLEDRefLowGainDiff[i])->SetBinContent(bin, diff);
966 ((TProfile*)fPedestalLEDRefLowGainDiff[i])->SetBinEntries(bin, 1);
967 ratio = GetPedLEDRefProfileLowGain(i)->GetBinContent(bin) / fReference->GetPedLEDRefProfileLowGain(i)->GetBinContent(bin);
968 ((TProfile*)fPedestalLEDRefLowGainRatio[i])->SetBinContent(bin, ratio);
969 ((TProfile*)fPedestalLEDRefLowGainRatio[i])->SetBinEntries(bin, 1);
978 //_____________________________________________________________________
979 void AliCaloCalibPedestal::ComputeHotAndWarningTowers(const char * hotMapFile)
980 { // look for hot/noisy towers
981 ValidateProfiles(); // make sure histos/profiles exist
983 char name[512];//Quite a long temp buffer, just in case the filename includes a path
986 snprintf(name, 512, "%s.txt", hotMapFile);
987 fout = new ofstream(name);
988 if (!fout->is_open()) {
990 fout = 0;//Set the pointer to empty if the file was not opened
994 for(int i = 0; i < fModules; i++){
996 //first we compute the peak-pedestal distribution for each supermodule...
997 if( GetPeakHighGainHisto(i)->GetEntries() > 0 ) {
998 double min = GetPeakProfileHighGain(i)->GetBinContent(GetPeakProfileHighGain(i)->GetMinimumBin());
999 double max = GetPeakProfileHighGain(i)->GetBinContent(GetPeakProfileHighGain(i)->GetMaximumBin());
1000 TH1D *hPeakFit = new TH1D(Form("hFit_%d", i), Form("hFit_%d", i), (int)((max-min)*10), min-1, max+1);
1002 for (int j = 1; j <= fColumns; j++) {
1003 for (int k = 1; k <= fRows; k++) {
1004 hPeakFit->Fill(GetPeakProfileHighGain(i)->GetBinContent(j, k));
1008 //...and then we try to fit it
1009 double mean = hPeakFit->GetMean();
1010 double sigma = hPeakFit->GetRMS();
1012 hPeakFit->Fit("gaus", "OQ", "", mean - 3*sigma, mean + 3*sigma);
1013 mean = hPeakFit->GetFunction("gaus")->GetParameter(1);
1014 sigma = hPeakFit->GetFunction("gaus")->GetParameter(2);
1016 catch (const std::exception & e) {
1017 printf("AliCaloCalibPedestal: TH1D PeakFit exception %s", e.what());
1023 //Then we look for warm/hot towers
1024 TH2F * hPeak2D = GetPeakHighGainHisto(i);
1025 hPeak2D->GetYaxis()->SetRangeUser( fWarningThreshold, hPeak2D->GetYaxis()->GetBinUpEdge(hPeak2D->GetNbinsY()) );
1028 int warnCounter = 0;
1029 for (int j = 1; j <= fColumns; j++) {
1030 for (int k = 1; k <= fRows; k++) {
1031 //we start looking for warm/warning towers...
1032 // histogram x-axis index
1033 idx = k-1 + fRows*(j-1); // this is what is used in the Fill call
1034 hPeak2D->GetXaxis()->SetRangeUser(idx, idx);
1035 warnCounter = (int) hPeak2D->Integral();
1036 if(warnCounter > fNEvents * fWarningFraction) {
1037 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kWarning);
1038 /* printf("mod %d col %d row %d warnCounter %d - status %d\n",
1039 i, j-1, k-1, warnCounter, (int) (kWarning)); */
1041 //...then we look for hot ones (towers whose values are greater than mean + X*sigma)
1042 if(GetPeakProfileHighGain(i)->GetBinContent(j, k) > mean + fHotSigma*sigma ) {
1043 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kHot);
1044 /* printf("mod %d col %d row %d binc %d - status %d\n",
1045 i, j-1, k-1, (int)(GetPeakProfileHighGain(i)->GetBinContent(j, k)), (int) (kHot)); */
1048 //Write the status to the hot/warm map file, if the file is open.
1049 // module - column - row - status (1=dead, 2= warm/warning , 3 = hot, see .h file enum)
1050 if (fout && ((TH2D*)fDeadMap[i])->GetBinContent(j, k) > 1) {
1055 << ((TH2D*)fDeadMap[i])->GetBinContent(j, k) << endl; }
1065 //_____________________________________________________________________
1066 void AliCaloCalibPedestal::ComputeDeadTowers(const char * deadMapFile)
1068 ValidateProfiles(); // make sure histos/profiles exist
1069 //Computes the number of dead towers etc etc into memory, after this you can call the GetDead... -functions
1073 ofstream * fout = 0;
1074 ofstream * diff = 0;
1075 char name[512];//Quite a long temp buffer, just in case the filename includes a path
1078 snprintf(name, 512, "%s.txt", deadMapFile);
1079 fout = new ofstream(name);
1080 snprintf(name, 512, "%sdiff.txt", deadMapFile);
1081 diff = new ofstream(name);
1082 if (!fout->is_open()) {
1084 fout = 0;//Set the pointer to empty if the file was not opened
1086 if (!diff->is_open()) {
1088 diff = 0;//Set the pointer to empty if the file was not opened
1092 for (int i = 0; i < fModules; i++) {
1093 if (GetPeakProfileHighGain(i)->GetEntries() > 0) { //don't care about empty histos
1094 for (int j = 1; j <= fColumns; j++) {
1095 for (int k = 1; k <= fRows; k++) {
1097 if (GetPeakProfileHighGain(i)->GetBinContent(j, k) < fDeadThreshold) {//It's dead
1098 countTot++;//One more dead total
1104 << "0" << endl;//Write the status to the deadmap file, if the file is open.
1107 if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) >= fDeadThreshold) {
1108 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kRecentlyDeceased);
1109 countNew++;//This tower wasn't dead before!
1111 ( *diff) << i << " "
1115 << "0" << endl;//Write the status to the deadmap difference file, if the file is open.
1119 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kDead);//This has been dead before. Nothing new
1122 else { //It's ALIVE!!
1123 //Don't bother with writing the live ones.
1124 if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) < fDeadThreshold) {
1125 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kResurrected);
1126 countRes++; //This tower was dead before => it's a miracle! :P
1132 << "1" << endl;//Write the status to the deadmap difference file, if the file is open.
1136 ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kAlive);
1141 }//end for j/columns
1142 }//end if GetEntries >= 0
1151 fDeadTowers = countTot;
1152 fNewDeadTowers = countNew;
1153 fResurrectedTowers = countRes;
1156 //_____________________________________________________________________
1157 Bool_t AliCaloCalibPedestal::IsBadChannel(int imod, int icol, int irow) const
1159 // Check if status info histo exists
1160 if (!fDeadMap[imod]) {
1164 //Check if channel is dead or hot.
1165 Int_t status = (Int_t) ( ((TH2D*)fDeadMap[imod])->GetBinContent(icol,irow) );
1166 if(status == kAlive)
1173 //_____________________________________________________________________
1174 void AliCaloCalibPedestal::SetChannelStatus(int imod, int icol, int irow, int status)
1176 ValidateProfiles(); // make sure histos/profiles exist
1177 //Set status of channel dead, hot, alive ...
1178 ((TH2D*)fDeadMap[imod])->SetBinContent(icol, irow, status);