X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliCaloCalibSignal.cxx;h=51d867073a0b5f5219c9569101d4bbeeacb2ba4a;hb=068e2c00b797fa2ccbb3a4852358a34ac05d32c5;hp=e77223e15339067f6ac7ed5aafa7c061af3466b0;hpb=f4fc542ceb7a9a4177394a4ec79eab52992f67c3;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliCaloCalibSignal.cxx b/EMCAL/AliCaloCalibSignal.cxx index e77223e1533..51d867073a0 100644 --- a/EMCAL/AliCaloCalibSignal.cxx +++ b/EMCAL/AliCaloCalibSignal.cxx @@ -26,18 +26,19 @@ */ // fed an event: // fSignals->ProcessEvent(fCaloRawStream,fRawEventHeaderBase); -// asked to draw graphs or profiles: -// fSignals->GetGraphAmpVsTimeHighGain(module,column,row)->Draw("ap"); -// or -// fSignals->GetProfAmpVsTimeHighGain(module,column,row)->Draw(); +// get some info: +// fSignals->GetXXX..() // etc. //________________________________________________________________________ +#include +#include +#include +#include "TProfile.h" #include "TFile.h" #include "AliRawReader.h" -#include "AliRawEventHeaderBase.h" -#include "AliCaloRawStream.h" +#include "AliCaloRawStreamV3.h" //The include file #include "AliCaloCalibSignal.h" @@ -46,6 +47,13 @@ ClassImp(AliCaloCalibSignal) using namespace std; +// variables for TTree filling; not sure if they should be static or not +static int fChannelNum; // for regular towers +static int fRefNum; // for LED +static double fAmp; +static double fAvgAmp; +static double fRMS; + // ctor; initialize everything in order to avoid compiler warnings // put some reasonable defaults AliCaloCalibSignal::AliCaloCalibSignal(kDetType detectorType) : @@ -53,25 +61,33 @@ AliCaloCalibSignal::AliCaloCalibSignal(kDetType detectorType) : fDetType(kNone), fColumns(0), fRows(0), + fLEDRefs(0), fModules(0), fCaloString(), fMapping(NULL), fRunNumber(-1), fStartTime(0), - fAmpCut(50), - fReqFractionAboveAmpCutVal(0.8), + fAmpCut(40), // min. 40 ADC counts as default + fReqFractionAboveAmpCutVal(0.6), // 60% in a strip, per default fReqFractionAboveAmp(kTRUE), + fAmpCutLEDRef(100), // min. 100 ADC counts as default + fReqLEDRefAboveAmpCutVal(kTRUE), fHour(0), fLatestHour(0), fUseAverage(kTRUE), fSecInAverage(1800), fNEvents(0), - fNAcceptedEvents(0) + fNAcceptedEvents(0), + fTreeAmpVsTime(NULL), + fTreeAvgAmpVsTime(NULL), + fTreeLEDAmpVsTime(NULL), + fTreeLEDAvgAmpVsTime(NULL) { //Default constructor. First we set the detector-type related constants. if (detectorType == kPhos) { fColumns = fgkPhosCols; fRows = fgkPhosRows; + fLEDRefs = fgkPhosLEDRefs; fModules = fgkPhosModules; fCaloString = "PHOS"; } @@ -79,47 +95,38 @@ AliCaloCalibSignal::AliCaloCalibSignal(kDetType detectorType) : //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 = fgkEmCalCols; - fRows = fgkEmCalRows; - fModules = fgkEmCalModules; + fColumns = AliEMCALGeoParams::fgkEMCALCols; + fRows = AliEMCALGeoParams::fgkEMCALRows; + fLEDRefs = AliEMCALGeoParams::fgkEMCALLEDRefs; + fModules = AliEMCALGeoParams::fgkEMCALModules; fCaloString = "EMCAL"; } fDetType = detectorType; - // Set the number of points for each Amp vs. Time graph to 0 - memset(fNHighGain, 0, sizeof(fNHighGain)); - memset(fNLowGain, 0, sizeof(fNLowGain)); - - CreateGraphs(); // set up the TGraphs - - // init TProfiles to NULL=0 also - memset(fProfAmpVsTimeHighGain, 0, sizeof(fProfAmpVsTimeHighGain)); - memset(fProfAmpVsTimeLowGain, 0, sizeof(fProfAmpVsTimeLowGain)); + ResetInfo(); // trees and counters } // dtor //_____________________________________________________________________ AliCaloCalibSignal::~AliCaloCalibSignal() { - ClearObjects(); + DeleteTrees(); } //_____________________________________________________________________ -void AliCaloCalibSignal::ClearObjects() +void AliCaloCalibSignal::DeleteTrees() { - // delete what was created in the ctor (TGraphs), and possible later (TProfiles) - for (int i=0; iSetName(name); - fGraphAmpVsTimeHighGain[id]->SetTitle(title); - fGraphAmpVsTimeHighGain[id]->SetMarkerStyle(20); - - // Low Gain - name = "fGraphAmpVsTimeLowGain_"; name += i; - name += "_"; name += ic; - name += "_"; name += ir; - title = "Amp vs. Time Low Gain Mod "; title += i; - title += " Col "; title += ic; - title += " Row "; title += ir; - - fGraphAmpVsTimeLowGain[id] = new TGraph(); - fGraphAmpVsTimeLowGain[id]->SetName(name); - fGraphAmpVsTimeLowGain[id]->SetTitle(title); - fGraphAmpVsTimeLowGain[id]->SetMarkerStyle(20); - - } - } + // initialize trees + // first, regular version + fTreeAmpVsTime = new TTree("fTreeAmpVsTime","Amplitude vs. Time Tree Variables"); + + fTreeAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I"); + fTreeAmpVsTime->Branch("fHour", &fHour, "fHour/D"); + fTreeAmpVsTime->Branch("fAmp", &fAmp, "fAmp/D"); + + // then, average version + fTreeAvgAmpVsTime = new TTree("fTreeAvgAmpVsTime","Average Amplitude vs. Time Tree Variables"); + + fTreeAvgAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I"); + fTreeAvgAmpVsTime->Branch("fHour", &fHour, "fHour/D"); + fTreeAvgAmpVsTime->Branch("fAvgAmp", &fAvgAmp, "fAvgAmp/D"); + fTreeAvgAmpVsTime->Branch("fRMS", &fRMS, "fRMS/D"); + + // then same for LED.. + fTreeLEDAmpVsTime = new TTree("fTreeLEDAmpVsTime","LED Amplitude vs. Time Tree Variables"); + fTreeLEDAmpVsTime->Branch("fRefNum", &fRefNum, "fRefNum/I"); + fTreeLEDAmpVsTime->Branch("fHour", &fHour, "fHour/D"); + fTreeLEDAmpVsTime->Branch("fAmp", &fAmp, "fAmp/D"); + + fTreeLEDAvgAmpVsTime = new TTree("fTreeLEDAvgAmpVsTime","Average LED Amplitude vs. Time Tree Variables"); + fTreeLEDAvgAmpVsTime->Branch("fRefNum", &fRefNum, "fRefNum/I"); + fTreeLEDAvgAmpVsTime->Branch("fHour", &fHour, "fHour/D"); + fTreeLEDAvgAmpVsTime->Branch("fAvgAmp", &fAvgAmp, "fAvgAmp/D"); + fTreeLEDAvgAmpVsTime->Branch("fRMS", &fRMS, "fRMS/D"); - }//end for nModules + return; } //_____________________________________________________________________ -void AliCaloCalibSignal::Reset() -{ +void AliCaloCalibSignal::ResetInfo() +{ // reset trees and counters Zero(); // set all counters to 0 - ClearObjects(); // delete previous TGraphs and TProfiles - CreateGraphs(); // and create some new ones + DeleteTrees(); // delete previous stuff + CreateTrees(); // and create some new ones return; } //_____________________________________________________________________ void AliCaloCalibSignal::Zero() { - // set all counters to 0; not cuts etc.though + // set all counters to 0; not cuts etc. though fHour = 0; fLatestHour = 0; fNEvents = 0; fNAcceptedEvents = 0; + + // Set the number of points for each tower: Amp vs. Time + memset(fNHighGain, 0, sizeof(fNHighGain)); + memset(fNLowGain, 0, sizeof(fNLowGain)); + // and LED reference + memset(fNRef, 0, sizeof(fNRef)); + return; } //_____________________________________________________________________ -Bool_t AliCaloCalibSignal::CheckFractionAboveAmp(int *AmpVal, int nTotChan) -{ - int nAbove = 0; +Bool_t AliCaloCalibSignal::CheckFractionAboveAmp(const int *iAmpVal, + int resultArray[]) +{ // check fraction of towers, per column, that are above amplitude cut + Bool_t returnCode = false; - int TowerNum = 0; + int iTowerNum = 0; + double fraction = 0; for (int i = 0; i fAmpCut) { + iTowerNum = GetTowerNum(i,j,k); + if (iAmpVal[iTowerNum] > fAmpCut) { nAbove++; } } + resultArray[i*fColumns +j] = 0; // init. to denied + if (nAbove > 0) { + fraction = (1.0*nAbove) / fRows; + /* + printf("DS mod %d col %d nAbove %d fraction %3.2f\n", + i, j, nAbove, fraction); + */ + if (fraction > fReqFractionAboveAmpCutVal) { + resultArray[i*fColumns + j] = nAbove; + returnCode = true; + } + } } - } + } // modules loop - double fraction = (1.0*nAbove) / nTotChan; + return returnCode; +} + + +//_____________________________________________________________________ +Bool_t AliCaloCalibSignal::CheckLEDRefAboveAmp(const int *iAmpVal, + int resultArray[]) +{ // check which LEDRef/Mon strips are above amplitude cut + Bool_t returnCode = false; + + int iRefNum = 0; + int gain = 1; // look at high gain; this should be rather saturated usually.. + for (int i = 0; i fAmpCutLEDRef) { + resultArray[i*fLEDRefs +j] = 1; // enough signal + returnCode = true; + } + else { + resultArray[i*fLEDRefs +j] = 0; // not enough signal + } + + /* + printf("DS mod %d LEDRef %d ampVal %d\n", + i, j, iAmpVal[iRefNum]); + */ + } // LEDRefs + } // modules loop - if (fraction > fReqFractionAboveAmpCutVal) { - return true; + return returnCode; +} + +// Parameter/cut handling +//_____________________________________________________________________ +void AliCaloCalibSignal::SetParametersFromFile(const char *parameterFile) +{ // set parameters from file + static const string delimitor("::"); + + // open, check input file + ifstream in( parameterFile ); + if( !in ) { + printf("in AliCaloCalibSignal::SetParametersFromFile - Using default/run_time parameters.\n"); + return; + } + + // 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). + + // 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 ) 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("AliCaloCalibSignal::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 + if ( (key == "fAmpCut") || (key == "fReqFractionAboveAmpCutVal") || + (key == "fAmpCutLEDRef") || (key == "fSecInAverage") ) { + istringstream iss(value); + printf("AliCaloCalibSignal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str()); + + if (key == "fAmpCut") { + iss >> fAmpCut; + } + else if (key == "fReqFractionAboveAmpCutVal") { + iss >> fReqFractionAboveAmpCutVal; + } + else if (key == "fAmpCutLEDRef") { + iss >> fAmpCutLEDRef; + } + else if (key == "fSecInAverage") { + iss >> fSecInAverage; + } + } // some match found/expected + + } } - else return false; + + in.close(); + return; +} + +//_____________________________________________________________________ +void AliCaloCalibSignal::WriteParametersToFile(const char *parameterFile) +{ // write parameters to file + static const string delimitor("::"); + ofstream out( parameterFile ); + out << "// " << parameterFile << endl; + out << "fAmpCut" << "::" << fAmpCut << endl; + out << "fReqFractionAboveAmpCutVal" << "::" << fReqFractionAboveAmpCutVal << endl; + out << "fAmpCutLEDRef" << "::" << fAmpCutLEDRef << endl; + out << "fSecInAverage" << "::" << fSecInAverage << endl; + + out.close(); + return; } //_____________________________________________________________________ Bool_t AliCaloCalibSignal::AddInfo(const AliCaloCalibSignal *sig) { - // just do this for the basic graphs/profiles that get filled in ProcessEvent - // may not have data for all channels, but let's just Add everything.. - // Note: this method will run into problems with TProfile adding if the binning of - // the local profiles is not the same as those provided by the argument *sig.. - int numGraphPoints = 0; - int id = 0; - int ip = 0; - for (int i = 0; i < fModules; i++) { - for (int j = 0; j < fColumns; j++) { - for (int k = 0; k < fRows; k++) { - - id = GetTowerNum(i,j,k); + // note/FIXME: we are not yet adding correctly the info for fN{HighGain,LowGain,Ref} here - but consider this a feature for now (20080905): we'll do Analyze() unless entries were found for a tower in this original object. + + // add info from sig's TTrees to ours.. + TTree *sigAmp = sig->GetTreeAmpVsTime(); + TTree *sigAvgAmp = sig->GetTreeAvgAmpVsTime(); + + // we could try some merging via TList or what also as a more elegant approach + // but I wanted with the stupid/simple and hopefully safe approach of looping + // over what we want to add.. + + // associate variables for sigAmp and sigAvgAmp: + sigAmp->SetBranchAddress("fChannelNum",&fChannelNum); + sigAmp->SetBranchAddress("fHour",&fHour); + sigAmp->SetBranchAddress("fAmp",&fAmp); + + // loop over the trees.. note that since we use the same variables we should not need + // to do any assignments between the getting and filling + for (int i=0; iGetEntries(); i++) { + sigAmp->GetEntry(i); + fTreeAmpVsTime->Fill(); + } - if(fUseAverage){ // add to Profiles - if (sig->GetProfAmpVsTimeHighGain(id)) { - GetProfAmpVsTimeHighGain(id)->Add(sig->GetProfAmpVsTimeHighGain(id)); - } - if (sig->GetProfAmpVsTimeLowGain(id)) { - GetProfAmpVsTimeLowGain(id)->Add(sig->GetProfAmpVsTimeLowGain(id)); - } - } - else{ // add to Graphs - // high gain - numGraphPoints= sig->GetGraphAmpVsTimeHighGain(id)->GetN(); - if (numGraphPoints > 0) { - // get the values - double *graphX = sig->GetGraphAmpVsTimeHighGain(id)->GetX(); - double *graphY = sig->GetGraphAmpVsTimeHighGain(id)->GetY(); - for(ip=0; ip < numGraphPoints; ip++){ - fGraphAmpVsTimeHighGain[id]->SetPoint(fNHighGain[id]++,graphX[ip],graphY[ip]); - } - } - // low gain - numGraphPoints= sig->GetGraphAmpVsTimeLowGain(id)->GetN(); - if (numGraphPoints > 0) { - // get the values - double *graphX = sig->GetGraphAmpVsTimeLowGain(id)->GetX(); - double *graphY = sig->GetGraphAmpVsTimeLowGain(id)->GetY(); - for(ip=0; ip < numGraphPoints; ip++){ - fGraphAmpVsTimeLowGain[id]->SetPoint(fNLowGain[id]++,graphX[ip],graphY[ip]); - } - } + sigAvgAmp->SetBranchAddress("fChannelNum",&fChannelNum); + sigAvgAmp->SetBranchAddress("fHour",&fHour); + sigAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp); + sigAvgAmp->SetBranchAddress("fRMS",&fRMS); - } + for (int i=0; iGetEntries(); i++) { + sigAvgAmp->GetEntry(i); + fTreeAvgAmpVsTime->Fill(); + } + + // also LED.. + TTree *sigLEDAmp = sig->GetTreeLEDAmpVsTime(); + TTree *sigLEDAvgAmp = sig->GetTreeLEDAvgAmpVsTime(); - }//end for nModules - }//end for nColumns - }//end for nRows + // associate variables for sigAmp and sigAvgAmp: + sigLEDAmp->SetBranchAddress("fRefNum",&fRefNum); + sigLEDAmp->SetBranchAddress("fHour",&fHour); + sigLEDAmp->SetBranchAddress("fAmp",&fAmp); + + // loop over the trees.. note that since we use the same variables we should not need + // to do any assignments between the getting and filling + for (int i=0; iGetEntries(); i++) { + sigLEDAmp->GetEntry(i); + fTreeLEDAmpVsTime->Fill(); + } - return kTRUE;//We succesfully added info from the supplied object + sigLEDAvgAmp->SetBranchAddress("fRefNum",&fRefNum); + sigLEDAvgAmp->SetBranchAddress("fHour",&fHour); + sigLEDAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp); + sigLEDAvgAmp->SetBranchAddress("fRMS",&fRMS); + + for (int i=0; iGetEntries(); i++) { + sigLEDAvgAmp->GetEntry(i); + fTreeLEDAvgAmpVsTime->Fill(); + } + + + return kTRUE;//We hopefully succesfully added info from the supplied object } //_____________________________________________________________________ Bool_t AliCaloCalibSignal::ProcessEvent(AliRawReader *rawReader) { // if fMapping is NULL the rawstream will crate its own mapping - AliCaloRawStream rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping); - - return ProcessEvent( &rawStream, (AliRawEventHeaderBase*)rawReader->GetEventHeader() ); + AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping); + if (fDetType == kEmCal) { + rawReader->Select("EMCAL", 0, AliEMCALGeoParams::fgkLastAltroDDL) ; //select EMCAL DDL range + } + return ProcessEvent( &rawStream, rawReader->GetTimestamp() ); } //_____________________________________________________________________ -Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStream *in, AliRawEventHeaderBase *aliHeader) +Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStreamV3 *in, UInt_t Timestamp) { // Method to process=analyze one event in the data stream if (!in) return kFALSE; //Return right away if there's a null pointer fNEvents++; // one more event - // PHOS has more towers than EMCAL, so use PHOS numbers to set array sizes - int AmpValHighGain[fgkMaxTowers]; - int AmpValLowGain[fgkMaxTowers]; + // use maximum numbers to set array sizes + int iAmpValHighGain[fgkMaxTowers]; + int iAmpValLowGain[fgkMaxTowers]; + memset(iAmpValHighGain, 0, sizeof(iAmpValHighGain)); + memset(iAmpValLowGain, 0, sizeof(iAmpValLowGain)); - memset(AmpValHighGain, 0, sizeof(AmpValHighGain)); - memset(AmpValLowGain, 0, sizeof(AmpValLowGain)); + // also for LED reference + int iLEDAmpVal[fgkMaxRefs * 2]; // factor 2 is for the two gain values + memset(iLEDAmpVal, 0, sizeof(iLEDAmpVal)); - int sample, i = 0; //The sample temp, and the sample number in current event. - int max = fgkSampleMin, min = fgkSampleMax;//Use these for picking the signal - int gain = 0; + int sample; // temporary value + int gain = 0; // high or low gain - // Number of Low and High gain channels for this event: + // Number of Low and High gain, and LED Ref, channels for this event: int nLowChan = 0; int nHighChan = 0; + int nLEDRefChan = 0; - int TowerNum = 0; // array index for TGraphs etc. + int iTowerNum = 0; // array index for regular towers + int iRefNum = 0; // array index for LED references // loop first to get the fraction of channels with amplitudes above cut - while (in->Next()) { - sample = in->GetSignal(); //Get the adc signal - if (sample < min) min = sample; - if (sample > max) max = sample; - i++; - if ( i >= in->GetTimeLength()) { + + while (in->NextDDL()) { + while (in->NextChannel()) { + + // counters + int max = AliEMCALGeoParams::fgkSampleMin, min = AliEMCALGeoParams::fgkSampleMax; // min and max sample values + int nsamples = 0; + + while (in->NextBunch()) { + const UShort_t *sig = in->GetSignals(); + nsamples += in->GetBunchLength(); + for (Int_t i = 0; i < in->GetBunchLength(); i++) { + sample = sig[i]; + + // check if it's a min or max value + if (sample < min) min = sample; + if (sample > max) max = sample; + + } // 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 + + gain = -1; // init to not valid value //If we're here then we're done with this tower - gain = 1 - in->IsLowGain(); - + if ( in->IsLowGain() ) { + gain = 0; + } + else if ( in->IsHighGain() ) { + gain = 1; + } + else if ( in->IsLEDMonData() ) { + gain = in->GetRow(); // gain coded in (in RCU/Altro mapping) as Row info for LED refs.. + } + + // 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); + printf("AliCaloCalibSignal::ProcessEvent = Oh no: arrayPos = %i.\n", arrayPos); + return kFALSE; } - // get tower number for AmpVal array - TowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow()); - - if (gain == 0) { - // fill amplitude into the array - AmpValLowGain[TowerNum] = max - min; - nLowChan++; - } - else if (gain==1) {//fill the high gain ones - // fill amplitude into the array - AmpValHighGain[TowerNum] = max - min; - nHighChan++; - }//end if gain - - - max = fgkSampleMin; min = fgkSampleMax; - i = 0; - - }//End if end of tower + if ( in->IsHighGain() || in->IsLowGain() ) { // regular tower + // get tower number for AmpVal array + iTowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow()); + + if (gain == 0) { + // fill amplitude into the array + iAmpValLowGain[iTowerNum] = max - min; + nLowChan++; + } + else if (gain==1) {//fill the high gain ones + // fill amplitude into the array + iAmpValHighGain[iTowerNum] = max - min; + nHighChan++; + }//end if gain + } // regular tower + else if ( in->IsLEDMonData() ) { // LED ref.; + // strip # is coded is 'column' in the channel maps + iRefNum = GetRefNum(arrayPos, in->GetColumn(), gain); + iLEDAmpVal[iRefNum] = max - min; + nLEDRefChan++; + } // end of LED ref + + } // nsamples>0 check, some data found for this channel; not only trailer/header + } // end while over channel - }//end while, of stream + }//end while over DDL's, of input stream - // now check if it was a led event, only use high gain (that should be sufficient) + in->Reset(); // just in case the next customer forgets to check if the stream was reset.. + + // now check if it was an LED event, using the LED Reference info per strip + + // by default all columns are accepted (init check to > 0) + int checkResultArray[AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALCols]; + for (int ia=0; ia<(AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALCols); ia++) { + checkResultArray[ia] = 1; + } if (fReqFractionAboveAmp) { bool ok = false; if (nHighChan > 0) { - ok = CheckFractionAboveAmp(AmpValHighGain, nHighChan); + ok = CheckFractionAboveAmp(iAmpValHighGain, checkResultArray); + } + if (!ok) return false; // skip event + } + + // by default all columns are accepted (init check to > 0) + int checkResultArrayLEDRef[AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALLEDRefs]; + for (int ia=0; ia<(AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALLEDRefs); ia++) { + checkResultArrayLEDRef[ia] = 1; + } + if (fReqLEDRefAboveAmpCutVal) { + bool ok = false; + if (nLEDRefChan > 0) { + ok = CheckLEDRefAboveAmp(iLEDAmpVal, checkResultArrayLEDRef); } if (!ok) return false; // skip event } @@ -395,68 +609,71 @@ Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStream *in, AliRawEventHeaderB fNAcceptedEvents++; // one more event accepted if (fStartTime == 0) { // if start-timestamp wasn't set,we'll pick it up from the first event we encounter - fStartTime = aliHeader->Get("Timestamp"); + fStartTime = Timestamp; } - fHour = (aliHeader->Get("Timestamp")-fStartTime)/(double)fgkNumSecInHr; + fHour = (Timestamp - fStartTime)/(double)fgkNumSecInHr; if (fLatestHour < fHour) { fLatestHour = fHour; } - // it is a led event, now fill graphs (maybe profiles later) - for(int i=0;i0 && checkResultArrayLEDRef[i*fLEDRefs + iStrip]>0) { // column passed check + for(int k=0; kSetPoint(fNHighGain[TowerNum]++,fHour,AmpValHighGain[TowerNum]); + if(iAmpValHighGain[iTowerNum]) { + fAmp = iAmpValHighGain[iTowerNum]; + fChannelNum = GetChannelNum(i,j,k,1); + fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValHighGain[iTowerNum]); + fNHighGain[iTowerNum]++; } - if(AmpValLowGain[TowerNum]) { - fGraphAmpVsTimeLowGain[TowerNum]->SetPoint(fNLowGain[TowerNum]++,fHour,AmpValLowGain[TowerNum]); + if(iAmpValLowGain[iTowerNum]) { + fAmp = iAmpValLowGain[iTowerNum]; + fChannelNum = GetChannelNum(i,j,k,0); + fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValLowGain[iTowerNum]); + fNLowGain[iTowerNum]++; } - } + } // rows + } // column passed check, and LED Ref for strip passed check (if any) + } // columns + + // also LED refs + for(int j=0; j0) || (checkResultArray[i*fColumns + iColFirst + 1]>0)) && // at least one column in strip passed check + (checkResultArrayLEDRef[i*fLEDRefs + j]>0) ) { // and LED Ref passed checks + for (gain=0; gain<2; gain++) { + fRefNum = GetRefNum(i, j, gain); + if (iLEDAmpVal[fRefNum]) { + fAmp = iLEDAmpVal[fRefNum]; + fTreeLEDAmpVsTime->Fill();//fRefNum,fHour,fAmp); + fNRef[fRefNum]++; + } + } // gain + } // at least one column in strip passed check, and LED Ref passed check (if any) } - } + + } // modules return kTRUE; } //_____________________________________________________________________ -void AliCaloCalibSignal::CreateProfile(int imod, int ic, int ir, int towerId, int gain, - int nbins, double min, double max) -{ //! create/setup a TProfile - TString title, name; - if (gain == 0) { - name = "fProfAmpVsTimeLowGain_"; - title = "Amp vs. Time Low Gain Mod "; - } - else if (gain == 1) { - name = "fProfAmpVsTimeHighGain_"; - title = "Amp vs. Time High Gain Mod "; - } - name += imod; - name += "_"; name += ic; - name += "_"; name += ir; - title += imod; - title += " Col "; title += ic; - title += " Row "; title += ir; - - // use "s" option for RMS - if (gain==0) { - fProfAmpVsTimeLowGain[towerId] = new TProfile(name,title, nbins, min, max,"s"); - } - else if (gain==1) { - fProfAmpVsTimeHighGain[towerId] = new TProfile(name,title, nbins, min, max,"s"); - } - - return; -} -//_____________________________________________________________________ -Bool_t AliCaloCalibSignal::Save(TString fileName, Bool_t saveEmptyGraphs) +Bool_t AliCaloCalibSignal::Save(TString fileName) { - //Saves all the histograms (or profiles, to be accurate) to the designated file + //Saves all the TTrees to the designated file TFile destFile(fileName, "recreate"); @@ -466,101 +683,187 @@ Bool_t AliCaloCalibSignal::Save(TString fileName, Bool_t saveEmptyGraphs) destFile.cd(); - // setup variables for the TProfile plot + // save the trees + fTreeAmpVsTime->Write(); + fTreeLEDAmpVsTime->Write(); + if (fUseAverage) { + Analyze(); // get the latest and greatest averages + fTreeAvgAmpVsTime->Write(); + fTreeLEDAvgAmpVsTime->Write(); + } + + destFile.Close(); + + return kTRUE; +} + +//_____________________________________________________________________ +Bool_t AliCaloCalibSignal::Analyze() +{ + // Fill the tree holding the average values + if (!fUseAverage) { return kFALSE; } + + // Reset the average TTree if Analyze has already been called earlier, + // meaning that the TTree could have been partially filled + if (fTreeAvgAmpVsTime->GetEntries() > 0) { + fTreeAvgAmpVsTime->Reset(); + } + + //0: setup variables for the TProfile plots that we'll use to do the averages int numProfBins = 0; double timeMin = 0; double timeMax = 0; - if (fUseAverage) { - if (fSecInAverage > 0) { - numProfBins = (int)( (fLatestHour*fgkNumSecInHr)/fSecInAverage + 1 ); // round-off - } - numProfBins += 2; // add extra buffer : first and last - double binSize = 1.0*fSecInAverage / fgkNumSecInHr; - timeMin = - binSize; - timeMax = timeMin + numProfBins*binSize; + if (fSecInAverage > 0) { + numProfBins = (int)( (fLatestHour*fgkNumSecInHr)/fSecInAverage + 1 ); // round-off } + numProfBins += 2; // add extra buffer : first and last + double binSize = 1.0*fSecInAverage / fgkNumSecInHr; + timeMin = - binSize; + timeMax = timeMin + numProfBins*binSize; - int numGraphPoints= 0; - int TowerNum = 0; - for (int i = 0; i < fModules; i++) { - - for(int ic=0;ic < fColumns;ic++){ - for(int ir=0;ir < fRows;ir++){ - - TowerNum = GetTowerNum(i, ic, ir); - - // 1st: high gain - numGraphPoints= fGraphAmpVsTimeHighGain[TowerNum]->GetN(); - if( numGraphPoints>0 || saveEmptyGraphs) { - - // average the graphs points over time if requested and put them in a profile plot - if(fUseAverage && numGraphPoints>0) { - - // get the values - double *graphX = fGraphAmpVsTimeHighGain[TowerNum]->GetX(); - double *graphY = fGraphAmpVsTimeHighGain[TowerNum]->GetY(); - - // create the TProfile: 1 is for High gain - CreateProfile(i, ic, ir, TowerNum, 1, - numProfBins, timeMin, timeMax); - - // loop over graph points and fill profile - for(int ip=0; ip < numGraphPoints; ip++){ - fProfAmpVsTimeHighGain[TowerNum]->Fill(graphX[ip],graphY[ip]); - } - - fProfAmpVsTimeHighGain[TowerNum]->GetXaxis()->SetTitle("Hours"); - fProfAmpVsTimeHighGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal"); - fProfAmpVsTimeHighGain[TowerNum]->Write(); + //1: set up TProfiles for the towers that had data + TProfile * profile[fgkMaxTowers*2]; // *2 is since we include both high and low gains + memset(profile, 0, sizeof(profile)); - } - else{ - //otherwise, just save the graphs and forget the profiling - fGraphAmpVsTimeHighGain[TowerNum]->GetXaxis()->SetTitle("Hours"); - fGraphAmpVsTimeHighGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal"); - fGraphAmpVsTimeHighGain[TowerNum]->Write(); - } - - } // low gain graph info should be saved in one form or another - - // 2nd: now go to the low gain case - numGraphPoints= fGraphAmpVsTimeLowGain[TowerNum]->GetN(); - if( numGraphPoints>0 || saveEmptyGraphs) { - - // average the graphs points over time if requested and put them in a profile plot - if(fUseAverage && numGraphPoints>0) { - - double *graphX = fGraphAmpVsTimeLowGain[TowerNum]->GetX(); - double *graphY = fGraphAmpVsTimeLowGain[TowerNum]->GetY(); - - // create the TProfile: 0 is for Low gain - CreateProfile(i, ic, ir, TowerNum, 0, - numProfBins, timeMin, timeMax); - - // loop over graph points and fill profile - for(int ip=0; ip < numGraphPoints; ip++){ - fProfAmpVsTimeLowGain[TowerNum]->Fill(graphX[ip],graphY[ip]); - } - - fProfAmpVsTimeLowGain[TowerNum]->GetXaxis()->SetTitle("Hours"); - fProfAmpVsTimeLowGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal"); - fProfAmpVsTimeLowGain[TowerNum]->Write(); + char name[200]; // for profile id and title + int iTowerNum = 0; - } - else{ - //otherwise, just save the graphs and forget the profiling - fGraphAmpVsTimeLowGain[TowerNum]->GetXaxis()->SetTitle("Hours"); - fGraphAmpVsTimeLowGain[TowerNum]->GetYaxis()->SetTitle("MaxAmplitude - Pedestal"); - fGraphAmpVsTimeLowGain[TowerNum]->Write(); - } - - } // low gain graph info should be saved in one form or another + for (int i = 0; i 0) { + fChannelNum = GetChannelNum(i, ic, ir, 1); + sprintf(name, "profileChan%d", fChannelNum); + profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s"); + } + + // same for low gain + if (fNLowGain[iTowerNum] > 0) { + fChannelNum = GetChannelNum(i, ic, ir, 0); + sprintf(name, "profileChan%d", fChannelNum); + profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s"); + } - } // fRows - } // fColumns + } // rows + } // columns + } // modules + + //2: fill profiles by looping over tree + // Set addresses for tree-readback also + fTreeAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum); + fTreeAmpVsTime->SetBranchAddress("fHour", &fHour); + fTreeAmpVsTime->SetBranchAddress("fAmp", &fAmp); + + for (int ient=0; ientGetEntries(); ient++) { + fTreeAmpVsTime->GetEntry(ient); + if (profile[fChannelNum]) { + // profile should always have been created above, for active channels + profile[fChannelNum]->Fill(fHour, fAmp); + } + } + + // re-associating the branch addresses here seems to be needed for OK 'average' storage + fTreeAvgAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum); + fTreeAvgAmpVsTime->SetBranchAddress("fHour", &fHour); + fTreeAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp); + fTreeAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS); + + //3: fill avg tree by looping over the profiles + for (fChannelNum = 0; fChannelNum<(fgkMaxTowers*2); fChannelNum++) { + if (profile[fChannelNum]) { // profile was created + if (profile[fChannelNum]->GetEntries() > 0) { // profile had some entries + for(int it=0; itGetBinEntries(it+1) > 0) { + fAvgAmp = profile[fChannelNum]->GetBinContent(it+1); + fHour = profile[fChannelNum]->GetBinCenter(it+1); + fRMS = profile[fChannelNum]->GetBinError(it+1); + fTreeAvgAmpVsTime->Fill(); + } // some entries for this bin + } // loop over bins + } // some entries for this profile + } // profile exists + } // loop over all possible channels + + + // and finally, go through same exercise for LED also.. + + //1: set up TProfiles for the towers that had data + TProfile * profileLED[fgkMaxRefs*2]; // *2 is since we include both high and low gains + memset(profileLED, 0, sizeof(profileLED)); + + for (int i = 0; i 0) { + sprintf(name, "profileLEDRef%d", fRefNum); + profileLED[fRefNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s"); + } + }// gain + } + } // modules + + //2: fill profiles by looping over tree + // Set addresses for tree-readback also + fTreeLEDAmpVsTime->SetBranchAddress("fRefNum", &fRefNum); + fTreeLEDAmpVsTime->SetBranchAddress("fHour", &fHour); + fTreeLEDAmpVsTime->SetBranchAddress("fAmp", &fAmp); + + for (int ient=0; ientGetEntries(); ient++) { + fTreeLEDAmpVsTime->GetEntry(ient); + if (profileLED[fRefNum]) { + // profile should always have been created above, for active channels + profileLED[fRefNum]->Fill(fHour, fAmp); + } + } + + // re-associating the branch addresses here seems to be needed for OK 'average' storage + fTreeLEDAvgAmpVsTime->SetBranchAddress("fRefNum", &fRefNum); + fTreeLEDAvgAmpVsTime->SetBranchAddress("fHour", &fHour); + fTreeLEDAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp); + fTreeLEDAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS); + + //3: fill avg tree by looping over the profiles + for (fRefNum = 0; fRefNum<(fgkMaxRefs*2); fRefNum++) { + if (profileLED[fRefNum]) { // profile was created + if (profileLED[fRefNum]->GetEntries() > 0) { // profile had some entries + for(int it=0; itGetBinEntries(it+1) > 0) { + fAvgAmp = profileLED[fRefNum]->GetBinContent(it+1); + fHour = profileLED[fRefNum]->GetBinCenter(it+1); + fRMS = profileLED[fRefNum]->GetBinError(it+1); + fTreeLEDAvgAmpVsTime->Fill(); + } // some entries for this bin + } // loop over bins + } // some entries for this profile + } // profile exists + } // loop over all possible channels + + // OK, we're done.. - } // fModules - destFile.Close(); - return kTRUE; } + +//_____________________________________________________________________ +Bool_t AliCaloCalibSignal::DecodeChannelNum(const int chanId, + int *imod, int *icol, int *irow, int *igain) const +{ // return the module, column, row, and gain for a given channel number + *igain = chanId/(fModules*fColumns*fRows); + *imod = (chanId/(fColumns*fRows)) % fModules; + *icol = (chanId/fRows) % fColumns; + *irow = chanId % fRows; + return kTRUE; +} + +//_____________________________________________________________________ +Bool_t AliCaloCalibSignal::DecodeRefNum(const int refId, + int *imod, int *istripMod, int *igain) const +{ // return the module, stripModule, and gain for a given reference number + *igain = refId/(fModules*fLEDRefs); + *imod = (refId/(fLEDRefs)) % fModules; + *istripMod = refId % fLEDRefs; + return kTRUE; +}