/************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ //* $Id$ */ //________________________________________________________________________ // // A help class for monitoring and calibration tools: MOOD, AMORE etc., // It can be created and used a la (ctor): /* //Create the object for making the histograms fPedestals = new AliCaloCalibPedestal( fDetType ); // AliCaloCalibPedestal knows how many modules we have for PHOS or EMCAL fNumModules = fPedestals->GetModules(); */ // fed an event: // fPedestals->ProcessEvent(fCaloRawStream); // asked to draw histograms: // fPedestals->GetDeadMap(i)->Draw("col"); // or // fPedestals->GetPeakProfileHighGainRatio((i < fNumModules) ? i : fVisibleModule)->Draw("colz"); // etc. // The pseudo-code examples above were from the first implementation in MOOD (summer 2007). //________________________________________________________________________ //#include "TCanvas.h" #include "TH1.h" #include "TF1.h" #include "TFile.h" #include #include #include #include #include #include "AliRawReader.h" #include "AliCaloRawStreamV3.h" //The include file #include "AliCaloCalibPedestal.h" ClassImp(AliCaloCalibPedestal) using namespace std; // ctor; initialize everything in order to avoid compiler warnings AliCaloCalibPedestal::AliCaloCalibPedestal(kDetType detectorType) : TObject(), fPedestalLowGain(), fPedestalHighGain(), fPedestalLEDRefLowGain(), fPedestalLEDRefHighGain(), fPeakMinusPedLowGain(), fPeakMinusPedHighGain(), fPeakMinusPedHighGainHisto(), fPedestalLowGainDiff(), fPedestalHighGainDiff(), fPedestalLEDRefLowGainDiff(), fPedestalLEDRefHighGainDiff(), fPeakMinusPedLowGainDiff(), fPeakMinusPedHighGainDiff(), fPedestalLowGainRatio(), fPedestalHighGainRatio(), fPedestalLEDRefLowGainRatio(), fPedestalLEDRefHighGainRatio(), fPeakMinusPedLowGainRatio(), fPeakMinusPedHighGainRatio(), fDeadMap(), fNEvents(0), fNChanFills(0), fDeadTowers(0), fNewDeadTowers(0), fResurrectedTowers(0), fReference(0), fDetType(kNone), fColumns(0), fRows(0), fLEDRefs(0), fModules(0), fRowMin(0), fRowMax(0), fRowMultiplier(0), fCaloString(), fMapping(NULL), fRunNumber(-1), fSelectPedestalSamples(kTRUE), fFirstPedestalSample(0), fLastPedestalSample(15), fDeadThreshold(5), fWarningThreshold(50), fWarningFraction(0.002), fHotSigma(5) { //Default constructor. First we set the detector-type related constants. if (detectorType == kPhos) { fColumns = fgkPhosCols; fRows = fgkPhosRows; fLEDRefs = fgkPhosLEDRefs; fModules = fgkPhosModules; fCaloString = "PHOS"; fRowMin = -1*fRows; fRowMax = 0; fRowMultiplier = -1; } else { //We'll just trust the enum to keep everything in line, so that if detectorType //isn't kPhos then it is kEmCal. Note, however, that this is not necessarily the //case, if someone intentionally gives another number fColumns = AliEMCALGeoParams::fgkEMCALCols; fRows = AliEMCALGeoParams::fgkEMCALRows; fLEDRefs = AliEMCALGeoParams::fgkEMCALLEDRefs; fModules = AliEMCALGeoParams::fgkEMCALModules; fCaloString = "EMCAL"; fRowMin = 0; fRowMax = fRows; fRowMultiplier = 1; } fDetType = detectorType; // ValidateProfiles(); // not to be done in ctor; info from Axel N. } //_____________________________________________________________________ void AliCaloCalibPedestal::ValidateProfiles() { //Make sure the basic histos exist if (!fPedestalLowGain.IsEmpty()) return; //The profiles already exist. We just check one, because they're all created at //the same time //Then, loop for the requested number of modules TString title, name; for (int i = 0; i < fModules; i++) { //Pedestals, low gain name = "hPedlowgain"; name += i; title = "Pedestals, low gain, module "; title += i; fPedestalLowGain.Add(new TProfile2D(name, title, fColumns, 0.0, fColumns, fRows, fRowMin, fRowMax,"s")); //Pedestals, high gain name = "hPedhighgain"; name += i; title = "Pedestals, high gain, module "; title += i; fPedestalHighGain.Add(new TProfile2D(name, title, fColumns, 0.0, fColumns, fRows, fRowMin, fRowMax,"s")); //LED Ref/Mon pedestals, low gain name = "hPedestalLEDReflowgain"; name += i; title = "Pedestal LEDRef, low gain, module "; title += i; fPedestalLEDRefLowGain.Add(new TProfile(name, title, fLEDRefs, 0.0, fLEDRefs, "s")); //LED Ref/Mon pedestals, high gain name = "hPedestalLEDRefhighgain"; name += i; title = "Pedestal LEDRef, high gain, module "; title += i; fPedestalLEDRefHighGain.Add(new TProfile(name, title, fLEDRefs, 0.0, fLEDRefs, "s")); //Peak-Pedestals, low gain name = "hPeakMinusPedlowgain"; name += i; title = "Peak-Pedestal, low gain, module "; title += i; fPeakMinusPedLowGain.Add(new TProfile2D(name, title, fColumns, 0.0, fColumns, fRows, fRowMin, fRowMax,"s")); //Peak-Pedestals, high gain name = "hPeakMinusPedhighgain"; name += i; title = "Peak-Pedestal, high gain, module "; title += i; fPeakMinusPedHighGain.Add(new TProfile2D(name, title, fColumns, 0.0, fColumns, fRows, fRowMin, fRowMax,"s")); //Peak-Pedestals, high gain - TH2F histo name = "hPeakMinusPedhighgainHisto"; name += i; title = "Peak-Pedestal, high gain, module "; title += i; fPeakMinusPedHighGainHisto.Add(new TH2F(name, title, fColumns*fRows, 0.0, fColumns*fRows, 100, 0, 1000)); name = "hDeadMap"; name += i; title = "Dead map, module "; title += i; fDeadMap.Add(new TH2D(name, title, fColumns, 0.0, fColumns, fRows, fRowMin, fRowMax)); }//end for nModules create the histograms CompressAndSetOwner(); } //_____________________________________________________________________ void AliCaloCalibPedestal::CompressAndSetOwner() { //Compress the arrays, in order to remove the empty objects (a 16 slot array is created by default) fPedestalLowGain.Compress(); fPedestalHighGain.Compress(); fPedestalLEDRefLowGain.Compress(); fPedestalLEDRefHighGain.Compress(); fPeakMinusPedLowGain.Compress(); fPeakMinusPedHighGain.Compress(); fPeakMinusPedHighGainHisto.Compress(); fDeadMap.Compress(); // set owner ship for everyone fPedestalLowGain.SetOwner(kTRUE); fPedestalHighGain.SetOwner(kTRUE); fPedestalLEDRefLowGain.SetOwner(kTRUE); fPedestalLEDRefHighGain.SetOwner(kTRUE); fPeakMinusPedLowGain.SetOwner(kTRUE); fPeakMinusPedHighGain.SetOwner(kTRUE); fPeakMinusPedHighGainHisto.SetOwner(kTRUE); fPedestalLowGainDiff.SetOwner(kTRUE); fPedestalHighGainDiff.SetOwner(kTRUE); fPedestalLEDRefLowGainDiff.SetOwner(kTRUE); fPedestalLEDRefHighGainDiff.SetOwner(kTRUE); fPeakMinusPedLowGainDiff.SetOwner(kTRUE); fPeakMinusPedHighGainDiff.SetOwner(kTRUE); fPedestalLowGainRatio.SetOwner(kTRUE); fPedestalHighGainRatio.SetOwner(kTRUE); fPedestalLEDRefLowGainRatio.SetOwner(kTRUE); fPedestalLEDRefHighGainRatio.SetOwner(kTRUE); fPeakMinusPedLowGainRatio.SetOwner(kTRUE); fPeakMinusPedHighGainRatio.SetOwner(kTRUE); fDeadMap.SetOwner(kTRUE); } // dtor //_____________________________________________________________________ AliCaloCalibPedestal::~AliCaloCalibPedestal() { //dtor if (fReference) delete fReference;//Delete the reference object, if it has been loaded // delete also TObjArray's fPedestalLowGain.Delete(); fPedestalHighGain.Delete(); fPedestalLEDRefLowGain.Delete(); fPedestalLEDRefHighGain.Delete(); fPeakMinusPedLowGain.Delete(); fPeakMinusPedHighGain.Delete(); fPeakMinusPedHighGainHisto.Delete(); fPedestalLowGainDiff.Delete(); fPedestalHighGainDiff.Delete(); fPedestalLEDRefLowGainDiff.Delete(); fPedestalLEDRefHighGainDiff.Delete(); fPeakMinusPedLowGainDiff.Delete(); fPeakMinusPedHighGainDiff.Delete(); fPedestalLowGainRatio.Delete(); fPedestalHighGainRatio.Delete(); fPedestalLEDRefLowGainRatio.Delete(); fPedestalLEDRefHighGainRatio.Delete(); fPeakMinusPedLowGainRatio.Delete(); fPeakMinusPedHighGainRatio.Delete(); fDeadMap.Delete(); } // copy ctor //_____________________________________________________________________ AliCaloCalibPedestal::AliCaloCalibPedestal(AliCaloCalibPedestal &ped) : TObject(ped), fPedestalLowGain(), fPedestalHighGain(), fPedestalLEDRefLowGain(), fPedestalLEDRefHighGain(), fPeakMinusPedLowGain(), fPeakMinusPedHighGain(), fPeakMinusPedHighGainHisto(), fPedestalLowGainDiff(), fPedestalHighGainDiff(), fPedestalLEDRefLowGainDiff(), fPedestalLEDRefHighGainDiff(), fPeakMinusPedLowGainDiff(), fPeakMinusPedHighGainDiff(), fPedestalLowGainRatio(), fPedestalHighGainRatio(), fPedestalLEDRefLowGainRatio(), fPedestalLEDRefHighGainRatio(), fPeakMinusPedLowGainRatio(), fPeakMinusPedHighGainRatio(), fDeadMap(), fNEvents(ped.GetNEvents()), fNChanFills(ped.GetNChanFills()), fDeadTowers(ped.GetDeadTowerCount()), fNewDeadTowers(ped.GetDeadTowerNew()), fResurrectedTowers(ped.GetDeadTowerResurrected()), fReference( 0 ), //! note that we do not try to copy the reference info here fDetType(ped.GetDetectorType()), fColumns(ped.GetColumns()), fRows(ped.GetRows()), fLEDRefs(ped.GetLEDRefs()), fModules(ped.GetModules()), fRowMin(ped.GetRowMin()), fRowMax(ped.GetRowMax()), fRowMultiplier(ped.GetRowMultiplier()), fCaloString(ped.GetCaloString()), fMapping(NULL), //! note that we are not copying the map info fRunNumber(ped.GetRunNumber()), fSelectPedestalSamples(ped.GetSelectPedestalSamples()), fFirstPedestalSample(ped.GetFirstPedestalSample()), fLastPedestalSample(ped.GetLastPedestalSample()), fDeadThreshold(ped.GetDeadThreshold()), fWarningThreshold(ped.GetWarningThreshold()), fWarningFraction(ped.GetWarningFraction()), fHotSigma(ped.GetHotSigma()) { // Then the ObjArray ones; we add the histograms rather than trying TObjArray = assignment //DS: this has not really been tested yet.. for (int i = 0; i < fModules; i++) { fPedestalLowGain.Add( ped.GetPedProfileLowGain(i) ); fPedestalHighGain.Add( ped.GetPedProfileHighGain(i) ); fPedestalLEDRefLowGain.Add( ped.GetPedLEDRefProfileLowGain(i) ); fPedestalLEDRefHighGain.Add( ped.GetPedLEDRefProfileHighGain(i) ); fPeakMinusPedLowGain.Add( ped.GetPeakProfileLowGain(i) ); fPeakMinusPedHighGain.Add( ped.GetPeakProfileHighGain(i) ); fPeakMinusPedHighGainHisto.Add( ped.GetPeakHighGainHisto(i) ); fDeadMap.Add( ped.GetDeadMap(i) ); }//end for nModules CompressAndSetOwner(); } // assignment operator; use copy ctor to make life easy.. //_____________________________________________________________________ AliCaloCalibPedestal& AliCaloCalibPedestal::operator = (AliCaloCalibPedestal &source) { // assignment operator; use copy ctor if (&source == this) return *this; new (this) AliCaloCalibPedestal(source); return *this; } //_____________________________________________________________________ void AliCaloCalibPedestal::Reset() { // Reset all arrays/histograms ValidateProfiles(); // make sure histos/profiles exist for (int i = 0; i < fModules; i++) { GetPedProfileLowGain(i)->Reset(); GetPedProfileHighGain(i)->Reset(); GetPedLEDRefProfileLowGain(i)->Reset(); GetPedLEDRefProfileHighGain(i)->Reset(); GetPeakProfileLowGain(i)->Reset(); GetPeakProfileHighGain(i)->Reset(); GetPeakHighGainHisto(i)->Reset(); GetDeadMap(i)->Reset(); if (!fPedestalLowGainDiff.IsEmpty()) { //This means that the comparison profiles have been created. GetPedProfileLowGainDiff(i)->Reset(); GetPedProfileHighGainDiff(i)->Reset(); GetPedLEDRefProfileLowGainDiff(i)->Reset(); GetPedLEDRefProfileHighGainDiff(i)->Reset(); GetPeakProfileLowGainDiff(i)->Reset(); GetPeakProfileHighGainDiff(i)->Reset(); GetPedProfileLowGainRatio(i)->Reset(); GetPedProfileHighGainRatio(i)->Reset(); GetPedLEDRefProfileLowGainRatio(i)->Reset(); GetPedLEDRefProfileHighGainRatio(i)->Reset(); GetPeakProfileLowGainRatio(i)->Reset(); GetPeakProfileHighGainRatio(i)->Reset(); } } fNEvents = 0; fNChanFills = 0; fDeadTowers = 0; fNewDeadTowers = 0; fResurrectedTowers = 0; //To think about: should fReference be deleted too?... let's not do it this time, at least... } // Parameter/cut handling //_____________________________________________________________________ void AliCaloCalibPedestal::SetParametersFromFile(const char *parameterFile) { // Note: this method is a bit more complicated than it really has to be // - allowing for multiple entries per line, arbitrary order of the // different variables etc. But I wanted to try and do this in as // correct a C++ way as I could (as an exercise). static const string delimitor("::"); // open, check input file ifstream in( parameterFile ); if( !in ) { printf("in AliCaloCalibPedestal::SetParametersFromFile - Using default/run_time parameters.\n"); return; } // read in char readline[1024]; while ((in.rdstate() & ios::failbit) == 0 ) { // Read into the raw char array and then construct a string // to do the searching in.getline(readline, 1024); istringstream s(readline); while ( ( s.rdstate() & ios::failbit ) == 0 ) { string keyValue; s >> keyValue; // check stream status if( ( s.rdstate() & ios::failbit ) == ios::failbit) break; // skip rest of line if comments found if( keyValue.substr( 0, 2 ) == "//" ) break; // look for "::" in keyValue pair size_t position = keyValue.find( delimitor ); if( position == string::npos ) { printf("wrong format for key::value pair: %s\n", keyValue.c_str()); } // split keyValue pair string key( keyValue.substr( 0, position ) ); string value( keyValue.substr( position+delimitor.size(), keyValue.size()-delimitor.size() ) ); // check value does not contain a new delimitor if( value.find( delimitor ) != string::npos ) { printf("wrong format for key::value pair: %s\n", keyValue.c_str()); } // debug: check key value pair // printf("AliCaloCalibPedestal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str()); // if the key matches with something we expect, we assign the new value istringstream iss(value); // the comparison strings defined at the beginning of this method if ( (key == "fFirstPedestalSample") || (key == "fLastPedestalSample") || (key == "fDeadThreshold") || (key == "fWarningThreshold") || (key == "fWarningFraction") || (key == "fHotSigma") ) { printf("AliCaloCalibPedestal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str()); if (key == "fFirstPedestalSample") { iss >> fFirstPedestalSample; } else if (key == "fLastPedestalSample") { iss >> fLastPedestalSample; } else if (key == "fDeadThreshold") { iss >> fDeadThreshold; } else if (key == "fWarningThreshold") { iss >> fWarningThreshold; } else if (key == "fWarningFraction") { iss >> fWarningFraction; } else if (key == "fHotSigma") { iss >> fHotSigma; } } // some match } } in.close(); return; } //_____________________________________________________________________ void AliCaloCalibPedestal::WriteParametersToFile(const char *parameterFile) { //Write parameters in file. static const string delimitor("::"); ofstream out( parameterFile ); out << "// " << parameterFile << endl; out << "fFirstPedestalSample" << "::" << fFirstPedestalSample << endl; out << "fLastPedestalSample" << "::" << fLastPedestalSample << endl; out << "fDeadThreshold" << "::" << fDeadThreshold << endl; out << "fWarningThreshold" << "::" << fWarningThreshold << endl; out << "fWarningFraction" << "::" << fWarningFraction << endl; out << "fHotSigma" << "::" << fHotSigma << endl; out.close(); return; } //_____________________________________________________________________ Bool_t AliCaloCalibPedestal::AddInfo(AliCaloCalibPedestal *ped) { // just do this for the basic histograms/profiles that get filled in ProcessEvent // may not have data for all modules, but let's just Add everything.. ValidateProfiles(); // make sure histos/profiles exist for (int i = 0; i < fModules; i++) { GetPedProfileLowGain(i)->Add( ped->GetPedProfileLowGain(i) ); GetPedProfileHighGain(i)->Add( ped->GetPedProfileHighGain(i) ); GetPedLEDRefProfileLowGain(i)->Add( ped->GetPedLEDRefProfileLowGain(i) ); GetPedLEDRefProfileHighGain(i)->Add( ped->GetPedLEDRefProfileHighGain(i) ); GetPeakProfileLowGain(i)->Add( ped->GetPeakProfileLowGain(i) ); GetPeakProfileHighGain(i)->Add( ped->GetPeakProfileHighGain(i) ); GetPeakHighGainHisto(i)->Add( ped->GetPeakHighGainHisto(i) ); }//end for nModules // We should also copy other pieces of info: counters and parameters // (not number of columns and rows etc which should be the same) // note that I just assign them here rather than Add them, but we // normally just Add (e.g. in Preprocessor) one object so this should be fine. fNEvents = ped->GetNEvents(); fNChanFills = ped->GetNChanFills(); fDeadTowers = ped->GetDeadTowerCount(); fNewDeadTowers = ped->GetDeadTowerNew(); fResurrectedTowers = ped->GetDeadTowerResurrected(); fRunNumber = ped->GetRunNumber(); fSelectPedestalSamples = ped->GetSelectPedestalSamples(); fFirstPedestalSample = ped->GetFirstPedestalSample(); fLastPedestalSample = ped->GetLastPedestalSample(); fDeadThreshold = ped->GetDeadThreshold(); fWarningThreshold = ped->GetWarningThreshold(); fWarningFraction = ped->GetWarningFraction(); fHotSigma = ped->GetHotSigma(); // DeadMap; Diff profiles etc would need to be redone after this operation return kTRUE;//We succesfully added info from the supplied object } //_____________________________________________________________________ Bool_t AliCaloCalibPedestal::ProcessEvent(AliRawReader *rawReader) { // if fMapping is NULL the rawstream will crate its own mapping AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping); if (fDetType == kEmCal) { rawReader->Select("EMCAL", 0, AliEMCALGeoParams::fgkLastAltroDDL) ; //select EMCAL DDL range } return ProcessEvent(&rawStream); } //_____________________________________________________________________ Bool_t AliCaloCalibPedestal::ProcessEvent(AliCaloRawStreamV3 *in) { // Method to process=analyze one event in the data stream if (!in) return kFALSE; //Return right away if there's a null pointer in->Reset(); // just in case the next customer forgets to check if the stream was reset.. fNEvents++; // one more event if (fNEvents==1) ValidateProfiles(); // 1st event, make sure histos/profiles exist // indices for the reading int sample = 0; int time = 0; int i = 0; // sample counter int startBin = 0; // start loop over input stream while (in->NextDDL()) { while (in->NextChannel()) { // counters int max = AliEMCALGeoParams::fgkSampleMin, min = AliEMCALGeoParams::fgkSampleMax; // min and max sample values int nsamples = 0; // pedestal samples int nPed = 0; vector pedSamples; while (in->NextBunch()) { const UShort_t *sig = in->GetSignals(); startBin = in->GetStartTimeBin(); nsamples += in->GetBunchLength(); for (i = 0; i < in->GetBunchLength(); i++) { sample = sig[i]; time = startBin--; // check if it's a min or max value if (sample < min) min = sample; if (sample > max) max = sample; // should we add it for the pedestal calculation? if ( (fFirstPedestalSample<=time && time<=fLastPedestalSample) || // sample time in range !fSelectPedestalSamples ) { // or we don't restrict the sample range.. - then we'll take all pedSamples.push_back( sig[i] ); nPed++; } } // loop over samples in bunch } // loop over bunches if (nsamples > 0) { // this check is needed for when we have zero-supp. on, but not sparse readout // it should be enough to check the SuperModule info for each DDL really, but let's keep it here for now int arrayPos = in->GetModule(); //The modules are numbered starting from 0 if (arrayPos >= fModules) { //TODO: return an error message, if appopriate (perhaps if debug>0?) return kFALSE; } //Debug if (arrayPos < 0 || arrayPos >= fModules) { printf("Oh no: arrayPos = %i.\n", arrayPos); } fNChanFills++; // one more channel found, and profile to be filled //NOTE: coordinates are (column, row) for the profiles if ( in->IsLowGain() ) { //fill the low gain histograms ((TProfile2D*)fPeakMinusPedLowGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min); if (nPed>0) { // only fill pedestal info in case it could be calculated for ( i=0; iFill(in->GetColumn(), fRowMultiplier*in->GetRow(), pedSamples[i]); } } } else if ( in->IsHighGain() ) { //fill the high gain ones ((TProfile2D*)fPeakMinusPedHighGain[arrayPos])->Fill(in->GetColumn(), fRowMultiplier*in->GetRow(), max - min); if (nPed>0) { // only fill pedestal info in case it could be calculated for ( i=0; iFill(in->GetColumn(), fRowMultiplier*in->GetRow(), pedSamples[i]); } } // for warning checks int idx = in->GetRow() + fRows * in->GetColumn(); ((TH2F*)fPeakMinusPedHighGainHisto[arrayPos])->Fill(idx, max - min); } else if ( in->IsLEDMonData() ) { // for LED Mon data, the mapping class holds the gain info in the Row variable // and the Strip number in the Column.. int gain = in->GetRow(); int stripId = in->GetColumn(); if (nPed>0 && stripIdFill(stripId, pedSamples[i]); } } else { for ( i=0; iFill(stripId, pedSamples[i]); } } } } } // nsamples>0 check, some data found for this channel; not only trailer/header }// end while over channel }//end while over DDL's, of input stream return kTRUE; } //_____________________________________________________________________ Bool_t AliCaloCalibPedestal::SaveHistograms(TString fileName, Bool_t saveEmptyHistos) { //Saves all the histograms (or profiles, to be accurate) to the designated file ValidateProfiles(); // make sure histos/profiles exist TFile destFile(fileName, "recreate"); if (destFile.IsZombie()) { return kFALSE; } destFile.cd(); for (int i = 0; i < fModules; i++) { if( ((TProfile2D *)fPeakMinusPedLowGain[i])->GetEntries() || saveEmptyHistos) { fPeakMinusPedLowGain[i]->Write(); } if( ((TProfile2D *)fPeakMinusPedHighGain[i])->GetEntries() || saveEmptyHistos) { fPeakMinusPedHighGain[i]->Write(); } if( ((TProfile2D *)fPedestalLowGain[i])->GetEntries() || saveEmptyHistos) { fPedestalLowGain[i]->Write(); } if( ((TProfile2D *)fPedestalHighGain[i])->GetEntries() || saveEmptyHistos) { fPedestalHighGain[i]->Write(); } if( ((TProfile *)fPedestalLEDRefLowGain[i])->GetEntries() || saveEmptyHistos) { fPedestalLEDRefLowGain[i]->Write(); } if( ((TProfile *)fPedestalLEDRefHighGain[i])->GetEntries() || saveEmptyHistos) { fPedestalLEDRefHighGain[i]->Write(); } if( ((TH2F *)fPeakMinusPedHighGainHisto[i])->GetEntries() || saveEmptyHistos) { fPeakMinusPedHighGainHisto[i]->Write(); } } destFile.Close(); return kTRUE; } //_____________________________________________________________________ Bool_t AliCaloCalibPedestal::LoadReferenceCalib(TString fileName, TString objectName) { //Make sure that the histograms created when loading the object are not destroyed as the file object is destroyed TH1::AddDirectory(kFALSE); TFile *sourceFile = new TFile(fileName); if (sourceFile->IsZombie()) { return kFALSE;//We couldn't load the reference } if (fReference) delete fReference;//Delete the reference object, if it already exists fReference = 0; fReference = (AliCaloCalibPedestal*)sourceFile->Get(objectName); if (!fReference || !(fReference->InheritsFrom(AliCaloCalibPedestal::Class())) || (fReference->GetDetectorType() != fDetType)) { if (fReference) delete fReference;//Delete the object, in case we had an object of the wrong type fReference = 0; return kFALSE; } delete sourceFile; //Reset the histogram ownership behaviour. NOTE: a better workaround would be good, since this may accidentally set AddDirectory to true, even //if we are called by someone who has set it to false... TH1::AddDirectory(kTRUE); return kTRUE;//We succesfully loaded the object } //_____________________________________________________________________ Bool_t AliCaloCalibPedestal::SetReference(AliCaloCalibPedestal *ref) { // set reference object if (fReference) delete fReference;//Delete the reference object, if it already exists fReference = 0; fReference = ref; if (!fReference || (fReference->GetDetectorType() != fDetType)) { if (fReference) delete fReference;//Delete the object, in case we had an object of the wrong type fReference = 0; return kFALSE; } return kTRUE;//We succesfully loaded the object } //_____________________________________________________________________ void AliCaloCalibPedestal::ValidateComparisonProfiles() { //Make sure the comparison histos exist if (!fPedestalLowGainDiff.IsEmpty()) return; //The profiles already exist. We just check one, because they're all created at //the same time //Then, loop for the requested number of modules TString title, name; for (int i = 0; i < fModules; i++) { //Pedestals, low gain name = "hPedlowgainDiff"; name += i; title = "Pedestals difference, low gain, module "; title += i; fPedestalLowGainDiff.Add(new TProfile2D(name, title, fColumns, 0.0, fColumns, fRows, fRowMin, fRowMax,"s")); //Pedestals, high gain name = "hPedhighgainDiff"; name += i; title = "Pedestals difference, high gain, module "; title += i; fPedestalHighGainDiff.Add(new TProfile2D(name, title, fColumns, 0.0, fColumns, fRows, fRowMin, fRowMax,"s")); //LED Ref/Mon pedestals, low gain name = "hPedestalLEDReflowgainDiff"; name += i; title = "Pedestal difference LEDRef, low gain, module "; title += i; fPedestalLEDRefLowGainDiff.Add(new TProfile(name, title, fLEDRefs, 0.0, fLEDRefs, "s")); //LED Ref/Mon pedestals, high gain name = "hPedestalLEDRefhighgainDiff"; name += i; title = "Pedestal difference LEDRef, high gain, module "; title += i; fPedestalLEDRefHighGainDiff.Add(new TProfile(name, title, fLEDRefs, 0.0, fLEDRefs, "s")); //Peak-Pedestals, high gain name = "hPeakMinusPedhighgainDiff"; name += i; title = "Peak-Pedestal difference, high gain, module "; title += i; fPeakMinusPedHighGainDiff.Add(new TProfile2D(name, title, fColumns, 0.0, fColumns, fRows, fRowMin, fRowMax,"s")); //Peak-Pedestals, low gain name = "hPeakMinusPedlowgainDiff"; name += i; title = "Peak-Pedestal difference, low gain, module "; title += i; fPeakMinusPedLowGainDiff.Add(new TProfile2D(name, title, fColumns, 0.0, fColumns, fRows, fRowMin, fRowMax,"s")); //Pedestals, low gain name = "hPedlowgainRatio"; name += i; title = "Pedestals ratio, low gain, module "; title += i; fPedestalLowGainRatio.Add(new TProfile2D(name, title, fColumns, 0.0, fColumns, fRows, fRowMin, fRowMax,"s")); //Pedestals, high gain name = "hPedhighgainRatio"; name += i; title = "Pedestals ratio, high gain, module "; title += i; fPedestalHighGainRatio.Add(new TProfile2D(name, title, fColumns, 0.0, fColumns, fRows, fRowMin, fRowMax,"s")); //LED Ref/Mon pedestals, low gain name = "hPedestalLEDReflowgainRatio"; name += i; title = "Pedestal ratio LEDRef, low gain, module "; title += i; fPedestalLEDRefLowGainRatio.Add(new TProfile(name, title, fLEDRefs, 0.0, fLEDRefs, "s")); //LED Ref/Mon pedestals, high gain name = "hPedestalLEDRefhighgainRatio"; name += i; title = "Pedestal ratio LEDRef, high gain, module "; title += i; fPedestalLEDRefHighGainRatio.Add(new TProfile(name, title, fLEDRefs, 0.0, fLEDRefs, "s")); //Peak-Pedestals, low gain name = "hPeakMinusPedlowgainRatio"; name += i; title = "Peak-Pedestal ratio, low gain, module "; title += i; fPeakMinusPedLowGainRatio.Add(new TProfile2D(name, title, fColumns, 0.0, fColumns, fRows, fRowMin, fRowMax,"s")); //Peak-Pedestals, high gain name = "hPeakMinusPedhighgainRatio"; name += i; title = "Peak-Pedestal ratio, high gain, module "; title += i; fPeakMinusPedHighGainRatio.Add(new TProfile2D(name, title, fColumns, 0.0, fColumns, fRows, fRowMin, fRowMax,"s")); }//end for nModules create the histograms } //_____________________________________________________________________ void AliCaloCalibPedestal::ComputeDiffAndRatio() { // calculate differences and ratios relative to a reference ValidateProfiles(); // make sure histos/profiles exist ValidateComparisonProfiles();//Make sure the comparison histos exist if (!fReference) { return;//Return if the reference object isn't loaded } int bin = 0; double diff = 0; double ratio = 1; for (int i = 0; i < fModules; i++) { //For computing the difference, we cannot simply do TProfile2D->Add(), because that subtracts the sum of all entries, //which means that the mean of the new profile will not be the difference of the means. So do it by hand: for (int j = 0; j < fColumns; j++) { for (int k = 0; k < fRows; k++) { bin = ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->GetBin(j+1, k+1);//Note that we assume here that all histos have the same structure... if (fReference->GetPeakProfileHighGain(i)->GetBinContent(bin) > 0) { diff = GetPeakProfileHighGain(i)->GetBinContent(bin) - fReference->GetPeakProfileHighGain(i)->GetBinContent(bin); ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinContent(bin, diff); ((TProfile2D*)fPeakMinusPedHighGainDiff[i])->SetBinEntries(bin, 1); ratio = GetPeakProfileHighGain(i)->GetBinContent(bin) / fReference->GetPeakProfileHighGain(i)->GetBinContent(bin); ((TProfile2D*)fPeakMinusPedHighGainRatio[i])->SetBinContent(bin, ratio); ((TProfile2D*)fPeakMinusPedHighGainRatio[i])->SetBinEntries(bin, 1); } if (fReference->GetPeakProfileLowGain(i)->GetBinContent(bin) > 0) { diff = GetPeakProfileLowGain(i)->GetBinContent(bin) - fReference->GetPeakProfileLowGain(i)->GetBinContent(bin); ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinContent(bin, diff); ((TProfile2D*)fPeakMinusPedLowGainDiff[i])->SetBinEntries(bin, 1); ratio = GetPeakProfileLowGain(i)->GetBinContent(bin) / fReference->GetPeakProfileLowGain(i)->GetBinContent(bin); ((TProfile2D*)fPeakMinusPedLowGainRatio[i])->SetBinContent(bin, ratio); ((TProfile2D*)fPeakMinusPedLowGainRatio[i])->SetBinEntries(bin, 1); } if (fReference->GetPedProfileHighGain(i)->GetBinContent(bin) > 0) { diff = GetPedProfileHighGain(i)->GetBinContent(bin) - fReference->GetPedProfileHighGain(i)->GetBinContent(bin); ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinContent(bin, diff); ((TProfile2D*)fPedestalHighGainDiff[i])->SetBinEntries(bin, 1); ratio = GetPedProfileHighGain(i)->GetBinContent(bin) / fReference->GetPedProfileHighGain(i)->GetBinContent(bin); ((TProfile2D*)fPedestalHighGainRatio[i])->SetBinContent(bin, ratio); ((TProfile2D*)fPedestalHighGainRatio[i])->SetBinEntries(bin, 1); } if (fReference->GetPedProfileLowGain(i)->GetBinContent(bin) > 0) { diff = GetPedProfileLowGain(i)->GetBinContent(bin) - fReference->GetPedProfileLowGain(i)->GetBinContent(bin); ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinContent(bin, diff); ((TProfile2D*)fPedestalLowGainDiff[i])->SetBinEntries(bin, 1); ratio = GetPedProfileLowGain(i)->GetBinContent(bin) / fReference->GetPedProfileLowGain(i)->GetBinContent(bin); ((TProfile2D*)fPedestalLowGainRatio[i])->SetBinContent(bin, ratio); ((TProfile2D*)fPedestalLowGainRatio[i])->SetBinEntries(bin, 1); } } // rows } // columns // same for LED Ref/Mon channels for (int j = 0; j <= fLEDRefs; j++) { bin = j+1;//Note that we assume here that all histos have the same structure... if (fReference->GetPedLEDRefProfileHighGain(i)->GetBinContent(bin) > 0) { diff = GetPedLEDRefProfileHighGain(i)->GetBinContent(bin) - fReference->GetPedLEDRefProfileHighGain(i)->GetBinContent(bin); ((TProfile*)fPedestalLEDRefHighGainDiff[i])->SetBinContent(bin, diff); ((TProfile*)fPedestalLEDRefHighGainDiff[i])->SetBinEntries(bin, 1); ratio = GetPedLEDRefProfileHighGain(i)->GetBinContent(bin) / fReference->GetPedLEDRefProfileHighGain(i)->GetBinContent(bin); ((TProfile*)fPedestalLEDRefHighGainRatio[i])->SetBinContent(bin, ratio); ((TProfile*)fPedestalLEDRefHighGainRatio[i])->SetBinEntries(bin, 1); } if (fReference->GetPedLEDRefProfileLowGain(i)->GetBinContent(bin) > 0) { diff = GetPedLEDRefProfileLowGain(i)->GetBinContent(bin) - fReference->GetPedLEDRefProfileLowGain(i)->GetBinContent(bin); ((TProfile*)fPedestalLEDRefLowGainDiff[i])->SetBinContent(bin, diff); ((TProfile*)fPedestalLEDRefLowGainDiff[i])->SetBinEntries(bin, 1); ratio = GetPedLEDRefProfileLowGain(i)->GetBinContent(bin) / fReference->GetPedLEDRefProfileLowGain(i)->GetBinContent(bin); ((TProfile*)fPedestalLEDRefLowGainRatio[i])->SetBinContent(bin, ratio); ((TProfile*)fPedestalLEDRefLowGainRatio[i])->SetBinEntries(bin, 1); } } } // modules } //_____________________________________________________________________ void AliCaloCalibPedestal::ComputeHotAndWarningTowers(const char * hotMapFile) { // look for hot/noisy towers ValidateProfiles(); // make sure histos/profiles exist ofstream * fout = 0; char name[512];//Quite a long temp buffer, just in case the filename includes a path if (hotMapFile) { snprintf(name, 512, "%s.txt", hotMapFile); fout = new ofstream(name); if (!fout->is_open()) { delete fout; fout = 0;//Set the pointer to empty if the file was not opened } } for(int i = 0; i < fModules; i++){ //first we compute the peak-pedestal distribution for each supermodule... if( GetPeakHighGainHisto(i)->GetEntries() > 0 ) { double min = GetPeakProfileHighGain(i)->GetBinContent(GetPeakProfileHighGain(i)->GetMinimumBin()); double max = GetPeakProfileHighGain(i)->GetBinContent(GetPeakProfileHighGain(i)->GetMaximumBin()); TH1D *hPeakFit = new TH1D(Form("hFit_%d", i), Form("hFit_%d", i), (int)((max-min)*10), min-1, max+1); for (int j = 1; j <= fColumns; j++) { for (int k = 1; k <= fRows; k++) { hPeakFit->Fill(GetPeakProfileHighGain(i)->GetBinContent(j, k)); } } //...and then we try to fit it double mean = hPeakFit->GetMean(); double sigma = hPeakFit->GetRMS(); try { hPeakFit->Fit("gaus", "OQ", "", mean - 3*sigma, mean + 3*sigma); mean = hPeakFit->GetFunction("gaus")->GetParameter(1); sigma = hPeakFit->GetFunction("gaus")->GetParameter(2); } catch (const std::exception & e) { printf("AliCaloCalibPedestal: TH1D PeakFit exception %s", e.what()); } //hPeakFit->Draw(); delete hPeakFit; //Then we look for warm/hot towers TH2F * hPeak2D = GetPeakHighGainHisto(i); hPeak2D->GetYaxis()->SetRangeUser( fWarningThreshold, hPeak2D->GetYaxis()->GetBinUpEdge(hPeak2D->GetNbinsY()) ); int idx = 0 ; int warnCounter = 0; for (int j = 1; j <= fColumns; j++) { for (int k = 1; k <= fRows; k++) { //we start looking for warm/warning towers... // histogram x-axis index idx = k-1 + fRows*(j-1); // this is what is used in the Fill call hPeak2D->GetXaxis()->SetRangeUser(idx, idx); warnCounter = (int) hPeak2D->Integral(); if(warnCounter > fNEvents * fWarningFraction) { ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kWarning); /* printf("mod %d col %d row %d warnCounter %d - status %d\n", i, j-1, k-1, warnCounter, (int) (kWarning)); */ } //...then we look for hot ones (towers whose values are greater than mean + X*sigma) if(GetPeakProfileHighGain(i)->GetBinContent(j, k) > mean + fHotSigma*sigma ) { ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kHot); /* printf("mod %d col %d row %d binc %d - status %d\n", i, j-1, k-1, (int)(GetPeakProfileHighGain(i)->GetBinContent(j, k)), (int) (kHot)); */ } //Write the status to the hot/warm map file, if the file is open. // module - column - row - status (1=dead, 2= warm/warning , 3 = hot, see .h file enum) if (fout && ((TH2D*)fDeadMap[i])->GetBinContent(j, k) > 1) { (*fout) << i << " " << (j - 1) << " " << (k - 1) << " " << ((TH2D*)fDeadMap[i])->GetBinContent(j, k) << endl; } } } } } return; } //_____________________________________________________________________ void AliCaloCalibPedestal::ComputeDeadTowers(const char * deadMapFile) { ValidateProfiles(); // make sure histos/profiles exist //Computes the number of dead towers etc etc into memory, after this you can call the GetDead... -functions int countTot = 0; int countNew = 0; int countRes = 0; ofstream * fout = 0; ofstream * diff = 0; char name[512];//Quite a long temp buffer, just in case the filename includes a path if (deadMapFile) { snprintf(name, 512, "%s.txt", deadMapFile); fout = new ofstream(name); snprintf(name, 512, "%sdiff.txt", deadMapFile); diff = new ofstream(name); if (!fout->is_open()) { delete fout; fout = 0;//Set the pointer to empty if the file was not opened } if (!diff->is_open()) { delete diff; diff = 0;//Set the pointer to empty if the file was not opened } } for (int i = 0; i < fModules; i++) { if (GetPeakProfileHighGain(i)->GetEntries() > 0) { //don't care about empty histos for (int j = 1; j <= fColumns; j++) { for (int k = 1; k <= fRows; k++) { if (GetPeakProfileHighGain(i)->GetBinContent(j, k) < fDeadThreshold) {//It's dead countTot++;//One more dead total if (fout) { (*fout) << i << " " << (j - 1) << " " << (k - 1) << " " << "1" << " " << "0" << endl;//Write the status to the deadmap file, if the file is open. } if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) >= fDeadThreshold) { ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kRecentlyDeceased); countNew++;//This tower wasn't dead before! if (diff) { ( *diff) << i << " " << (j - 1) << " " << (k - 1) << " " << "1" << " " << "0" << endl;//Write the status to the deadmap difference file, if the file is open. } } else { ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kDead);//This has been dead before. Nothing new } } else { //It's ALIVE!! //Don't bother with writing the live ones. if (fReference && fReference->GetPeakProfileHighGain(i)->GetBinContent(j, k) < fDeadThreshold) { ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kResurrected); countRes++; //This tower was dead before => it's a miracle! :P if (diff) { (*diff) << i << " " << (j - 1) << " " << (k - 1) << " " << "1" << " " << "1" << endl;//Write the status to the deadmap difference file, if the file is open. } } else { ((TH2D*)fDeadMap[i])->SetBinContent(j, k, kAlive); } } }//end for k/rows }//end for j/columns }//end if GetEntries >= 0 }//end for modules if (fout) { fout->close(); delete fout; } fDeadTowers = countTot; fNewDeadTowers = countNew; fResurrectedTowers = countRes; } //_____________________________________________________________________ Bool_t AliCaloCalibPedestal::IsBadChannel(int imod, int icol, int irow) const { // Check if status info histo exists if (!fDeadMap[imod]) { return kFALSE; } //Check if channel is dead or hot. Int_t status = (Int_t) ( ((TH2D*)fDeadMap[imod])->GetBinContent(icol,irow) ); if(status == kAlive) return kFALSE; else return kTRUE; } //_____________________________________________________________________ void AliCaloCalibPedestal::SetChannelStatus(int imod, int icol, int irow, int status) { ValidateProfiles(); // make sure histos/profiles exist //Set status of channel dead, hot, alive ... ((TH2D*)fDeadMap[imod])->SetBinContent(icol, irow, status); }