1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15 /* $Id: AliCaloCalibSignal.cxx $ */
17 //________________________________________________________________________
19 // A help class for monitoring and calibration tools: MOOD, AMORE etc.,
20 // It can be created and used a la (ctor):
22 //Create the object for making the histograms
23 fSignals = new AliCaloCalibSignal( fDetType );
24 // AliCaloCalibSignal knows how many modules we have for PHOS or EMCAL
25 fNumModules = fSignals->GetModules();
28 // fSignals->ProcessEvent(fCaloRawStream,fRawEventHeaderBase);
30 // fSignals->GetXXX..()
32 //________________________________________________________________________
40 #include "AliRawReader.h"
41 #include "AliCaloRawStreamV3.h"
44 #include "AliCaloCalibSignal.h"
46 ClassImp(AliCaloCalibSignal)
50 // variables for TTree filling; not sure if they should be static or not
51 static int fChannelNum; // for regular towers
52 static int fRefNum; // for LED
54 static double fAvgAmp;
57 // ctor; initialize everything in order to avoid compiler warnings
58 // put some reasonable defaults
59 AliCaloCalibSignal::AliCaloCalibSignal(kDetType detectorType) :
70 fAmpCut(40), // min. 40 ADC counts as default
71 fReqFractionAboveAmpCutVal(0.6), // 60% in a strip, per default
72 fReqFractionAboveAmp(kTRUE),
80 fTreeAvgAmpVsTime(NULL),
81 fTreeLEDAmpVsTime(NULL),
82 fTreeLEDAvgAmpVsTime(NULL)
84 //Default constructor. First we set the detector-type related constants.
85 if (detectorType == kPhos) {
86 fColumns = fgkPhosCols;
88 fLEDRefs = fgkPhosLEDRefs;
89 fModules = fgkPhosModules;
93 //We'll just trust the enum to keep everything in line, so that if detectorType
94 //isn't kPhos then it is kEmCal. Note, however, that this is not necessarily the
95 //case, if someone intentionally gives another number
96 fColumns = AliEMCALGeoParams::fgkEMCALCols;
97 fRows = AliEMCALGeoParams::fgkEMCALRows;
98 fLEDRefs = AliEMCALGeoParams::fgkEMCALLEDRefs;
99 fModules = AliEMCALGeoParams::fgkEMCALModules;
100 fCaloString = "EMCAL";
103 fDetType = detectorType;
105 ResetInfo(); // trees and counters
109 //_____________________________________________________________________
110 AliCaloCalibSignal::~AliCaloCalibSignal()
115 //_____________________________________________________________________
116 void AliCaloCalibSignal::DeleteTrees()
118 // delete what was created in the ctor (TTrees)
119 if (fTreeAmpVsTime) delete fTreeAmpVsTime;
120 if (fTreeAvgAmpVsTime) delete fTreeAvgAmpVsTime;
121 if (fTreeLEDAmpVsTime) delete fTreeLEDAmpVsTime;
122 if (fTreeLEDAvgAmpVsTime) delete fTreeLEDAvgAmpVsTime;
123 // and reset pointers
124 fTreeAmpVsTime = NULL;
125 fTreeAvgAmpVsTime = NULL;
126 fTreeLEDAmpVsTime = NULL;
127 fTreeLEDAvgAmpVsTime = NULL;
133 //_____________________________________________________________________
134 AliCaloCalibSignal::AliCaloCalibSignal(const AliCaloCalibSignal &sig) :
136 fDetType(sig.GetDetectorType()),
137 fColumns(sig.GetColumns()),
138 fRows(sig.GetRows()),
139 fLEDRefs(sig.GetLEDRefs()),
140 fModules(sig.GetModules()),
141 fCaloString(sig.GetCaloString()),
142 fMapping(NULL), //! note that we are not copying the map info
143 fRunNumber(sig.GetRunNumber()),
144 fStartTime(sig.GetStartTime()),
145 fAmpCut(sig.GetAmpCut()),
146 fReqFractionAboveAmpCutVal(sig.GetReqFractionAboveAmpCutVal()),
147 fReqFractionAboveAmp(sig.GetReqFractionAboveAmp()),
148 fHour(sig.GetHour()),
149 fLatestHour(sig.GetLatestHour()),
150 fUseAverage(sig.GetUseAverage()),
151 fSecInAverage(sig.GetSecInAverage()),
152 fNEvents(sig.GetNEvents()),
153 fNAcceptedEvents(sig.GetNAcceptedEvents()),
154 fTreeAmpVsTime(NULL),
155 fTreeAvgAmpVsTime(NULL),
156 fTreeLEDAmpVsTime(NULL),
157 fTreeLEDAvgAmpVsTime(NULL)
159 // also the TTree contents
163 // assignment operator; use copy ctor to make life easy..
164 //_____________________________________________________________________
165 AliCaloCalibSignal& AliCaloCalibSignal::operator = (const AliCaloCalibSignal &source)
167 // assignment operator; use copy ctor
168 if (&source == this) return *this;
170 new (this) AliCaloCalibSignal(source);
174 //_____________________________________________________________________
175 void AliCaloCalibSignal::CreateTrees()
178 // first, regular version
179 fTreeAmpVsTime = new TTree("fTreeAmpVsTime","Amplitude vs. Time Tree Variables");
181 fTreeAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I");
182 fTreeAmpVsTime->Branch("fHour", &fHour, "fHour/D");
183 fTreeAmpVsTime->Branch("fAmp", &fAmp, "fAmp/D");
185 // then, average version
186 fTreeAvgAmpVsTime = new TTree("fTreeAvgAmpVsTime","Average Amplitude vs. Time Tree Variables");
188 fTreeAvgAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I");
189 fTreeAvgAmpVsTime->Branch("fHour", &fHour, "fHour/D");
190 fTreeAvgAmpVsTime->Branch("fAvgAmp", &fAvgAmp, "fAvgAmp/D");
191 fTreeAvgAmpVsTime->Branch("fRMS", &fRMS, "fRMS/D");
193 // then same for LED..
194 fTreeLEDAmpVsTime = new TTree("fTreeLEDAmpVsTime","LED Amplitude vs. Time Tree Variables");
195 fTreeLEDAmpVsTime->Branch("fRefNum", &fRefNum, "fRefNum/I");
196 fTreeLEDAmpVsTime->Branch("fHour", &fHour, "fHour/D");
197 fTreeLEDAmpVsTime->Branch("fAmp", &fAmp, "fAmp/D");
199 fTreeLEDAvgAmpVsTime = new TTree("fTreeLEDAvgAmpVsTime","Average LED Amplitude vs. Time Tree Variables");
200 fTreeLEDAvgAmpVsTime->Branch("fRefNum", &fRefNum, "fRefNum/I");
201 fTreeLEDAvgAmpVsTime->Branch("fHour", &fHour, "fHour/D");
202 fTreeLEDAvgAmpVsTime->Branch("fAvgAmp", &fAvgAmp, "fAvgAmp/D");
203 fTreeLEDAvgAmpVsTime->Branch("fRMS", &fRMS, "fRMS/D");
208 //_____________________________________________________________________
209 void AliCaloCalibSignal::ResetInfo()
211 Zero(); // set all counters to 0
212 DeleteTrees(); // delete previous stuff
213 CreateTrees(); // and create some new ones
217 //_____________________________________________________________________
218 void AliCaloCalibSignal::Zero()
220 // set all counters to 0; not cuts etc. though
224 fNAcceptedEvents = 0;
226 // Set the number of points for each tower: Amp vs. Time
227 memset(fNHighGain, 0, sizeof(fNHighGain));
228 memset(fNLowGain, 0, sizeof(fNLowGain));
230 memset(fNRef, 0, sizeof(fNRef));
235 //_____________________________________________________________________
236 Bool_t AliCaloCalibSignal::CheckFractionAboveAmp(int *AmpVal,
239 Bool_t returnCode = false;
243 for (int i = 0; i<fModules; i++) {
244 for (int j = 0; j<fColumns; j++) {
246 for (int k = 0; k<fRows; k++) {
247 TowerNum = GetTowerNum(i,j,k);
248 if (AmpVal[TowerNum] > fAmpCut) {
252 resultArray[i*fColumns +j] = 0; // init. to denied
254 fraction = (1.0*nAbove) / fRows;
256 printf("DS mod %d col %d nAbove %d fraction %3.2f\n",
257 i, j, nAbove, fraction);
259 if (fraction > fReqFractionAboveAmpCutVal) {
260 resultArray[i*fColumns + j] = nAbove;
270 // Parameter/cut handling
271 //_____________________________________________________________________
272 void AliCaloCalibSignal::SetParametersFromFile(const char *parameterFile)
274 static const string delimitor("::");
276 // open, check input file
277 ifstream in( parameterFile );
279 printf("in AliCaloCalibSignal::SetParametersFromFile - Using default/run_time parameters.\n");
283 // Note: this method is a bit more complicated than it really has to be
284 // - allowing for multiple entries per line, arbitrary order of the
285 // different variables etc. But I wanted to try and do this in as
286 // correct a C++ way as I could (as an exercise).
290 while ((in.rdstate() & ios::failbit) == 0 ) {
292 // Read into the raw char array and then construct a string
293 // to do the searching
294 in.getline(readline, 1024);
295 istringstream s(readline);
297 while ( ( s.rdstate() & ios::failbit ) == 0 ) {
302 // check stream status
303 if( s.rdstate() & ios::failbit ) break;
305 // skip rest of line if comments found
306 if( key_value.substr( 0, 2 ) == "//" ) break;
308 // look for "::" in key_value pair
309 size_t position = key_value.find( delimitor );
310 if( position == string::npos ) {
311 printf("wrong format for key::value pair: %s\n", key_value.c_str());
314 // split key_value pair
315 string key( key_value.substr( 0, position ) );
316 string value( key_value.substr( position+delimitor.size(),
317 key_value.size()-delimitor.size() ) );
319 // check value does not contain a new delimitor
320 if( value.find( delimitor ) != string::npos ) {
321 printf("wrong format for key::value pair: %s\n", key_value.c_str());
324 // debug: check key value pair
325 // printf("AliCaloCalibSignal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
327 // if the key matches with something we expect, we assign the new value
328 if ( (key == "fAmpCut") || (key == "fReqFractionAboveAmpCutVal") || (key == "fSecInAverage") ) {
329 istringstream iss(value);
330 printf("AliCaloCalibSignal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
332 if (key == "fAmpCut") {
335 else if (key == "fReqFractionAboveAmpCutVal") {
336 iss >> fReqFractionAboveAmpCutVal;
338 else if (key == "fSecInAverage") {
339 iss >> fSecInAverage;
341 } // some match found/expected
351 //_____________________________________________________________________
352 void AliCaloCalibSignal::WriteParametersToFile(const char *parameterFile)
354 static const string delimitor("::");
355 ofstream out( parameterFile );
356 out << "// " << parameterFile << endl;
357 out << "fAmpCut" << "::" << fAmpCut << endl;
358 out << "fReqFractionAboveAmpCutVal" << "::" << fReqFractionAboveAmpCutVal << endl;
359 out << "fSecInAverage" << "::" << fSecInAverage << endl;
365 //_____________________________________________________________________
366 Bool_t AliCaloCalibSignal::AddInfo(const AliCaloCalibSignal *sig)
368 // 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.
370 // add info from sig's TTrees to ours..
371 TTree *sigAmp = sig->GetTreeAmpVsTime();
372 TTree *sigAvgAmp = sig->GetTreeAvgAmpVsTime();
374 // we could try some merging via TList or what also as a more elegant approach
375 // but I wanted with the stupid/simple and hopefully safe approach of looping
376 // over what we want to add..
378 // associate variables for sigAmp and sigAvgAmp:
379 sigAmp->SetBranchAddress("fChannelNum",&fChannelNum);
380 sigAmp->SetBranchAddress("fHour",&fHour);
381 sigAmp->SetBranchAddress("fAmp",&fAmp);
383 // loop over the trees.. note that since we use the same variables we should not need
384 // to do any assignments between the getting and filling
385 for (int i=0; i<sigAmp->GetEntries(); i++) {
387 fTreeAmpVsTime->Fill();
390 sigAvgAmp->SetBranchAddress("fChannelNum",&fChannelNum);
391 sigAvgAmp->SetBranchAddress("fHour",&fHour);
392 sigAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
393 sigAvgAmp->SetBranchAddress("fRMS",&fRMS);
395 for (int i=0; i<sigAvgAmp->GetEntries(); i++) {
396 sigAvgAmp->GetEntry(i);
397 fTreeAvgAmpVsTime->Fill();
401 TTree *sigLEDAmp = sig->GetTreeLEDAmpVsTime();
402 TTree *sigLEDAvgAmp = sig->GetTreeLEDAvgAmpVsTime();
404 // associate variables for sigAmp and sigAvgAmp:
405 sigLEDAmp->SetBranchAddress("fRefNum",&fRefNum);
406 sigLEDAmp->SetBranchAddress("fHour",&fHour);
407 sigLEDAmp->SetBranchAddress("fAmp",&fAmp);
409 // loop over the trees.. note that since we use the same variables we should not need
410 // to do any assignments between the getting and filling
411 for (int i=0; i<sigLEDAmp->GetEntries(); i++) {
412 sigLEDAmp->GetEntry(i);
413 fTreeLEDAmpVsTime->Fill();
416 sigLEDAvgAmp->SetBranchAddress("fRefNum",&fRefNum);
417 sigLEDAvgAmp->SetBranchAddress("fHour",&fHour);
418 sigLEDAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
419 sigLEDAvgAmp->SetBranchAddress("fRMS",&fRMS);
421 for (int i=0; i<sigLEDAvgAmp->GetEntries(); i++) {
422 sigLEDAvgAmp->GetEntry(i);
423 fTreeLEDAvgAmpVsTime->Fill();
427 return kTRUE;//We hopefully succesfully added info from the supplied object
430 //_____________________________________________________________________
431 Bool_t AliCaloCalibSignal::ProcessEvent(AliRawReader *rawReader)
433 // if fMapping is NULL the rawstream will crate its own mapping
434 AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
436 return ProcessEvent( &rawStream, rawReader->GetTimestamp() );
439 //_____________________________________________________________________
440 Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStreamV3 *in, UInt_t Timestamp)
442 // Method to process=analyze one event in the data stream
443 if (!in) return kFALSE; //Return right away if there's a null pointer
445 fNEvents++; // one more event
447 // use maximum numbers to set array sizes
448 int AmpValHighGain[fgkMaxTowers];
449 int AmpValLowGain[fgkMaxTowers];
450 memset(AmpValHighGain, 0, sizeof(AmpValHighGain));
451 memset(AmpValLowGain, 0, sizeof(AmpValLowGain));
453 // also for LED reference
454 int LEDAmpVal[fgkMaxRefs * 2]; // factor 2 is for the two gain values
455 memset(LEDAmpVal, 0, sizeof(LEDAmpVal));
457 int sample; // temporary value
458 int gain = 0; // high or low gain
460 // Number of Low and High gain channels for this event:
464 int TowerNum = 0; // array index for regular towers
465 int RefNum = 0; // array index for LED references
467 // loop first to get the fraction of channels with amplitudes above cut
469 while (in->NextDDL()) {
470 while (in->NextChannel()) {
473 int max = AliEMCALGeoParams::fgkSampleMin, min = AliEMCALGeoParams::fgkSampleMax; // min and max sample values
476 while (in->NextBunch()) {
477 const UShort_t *sig = in->GetSignals();
478 nsamples += in->GetBunchLength();
479 for (Int_t i = 0; i < in->GetBunchLength(); i++) {
482 // check if it's a min or max value
483 if (sample < min) min = sample;
484 if (sample > max) max = sample;
486 } // loop over samples in bunch
487 } // loop over bunches
489 if (nsamples > 0) { // this check is needed for when we have zero-supp. on, but not sparse readout
491 gain = -1; // init to not valid value
492 //If we're here then we're done with this tower
493 if ( in->IsLowGain() ) {
496 else if ( in->IsHighGain() ) {
499 else if ( in->IsLEDMonData() ) {
500 gain = in->GetRow(); // gain coded in (in RCU/Altro mapping) as Row info for LED refs..
503 // it should be enough to check the SuperModule info for each DDL really, but let's keep it here for now
504 int arrayPos = in->GetModule(); //The modules are numbered starting from 0
506 if (arrayPos < 0 || arrayPos >= fModules) {
507 printf("AliCaloCalibSignal::ProcessEvent = Oh no: arrayPos = %i.\n", arrayPos);
511 if ( in->IsHighGain() || in->IsLowGain() ) { // regular tower
512 // get tower number for AmpVal array
513 TowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow());
516 // fill amplitude into the array
517 AmpValLowGain[TowerNum] = max - min;
520 else if (gain==1) {//fill the high gain ones
521 // fill amplitude into the array
522 AmpValHighGain[TowerNum] = max - min;
526 else if ( in->IsLEDMonData() ) { // LED ref.;
527 // strip # is coded is 'column' in the channel maps
528 RefNum = GetRefNum(arrayPos, in->GetColumn(), gain);
529 LEDAmpVal[RefNum] = max - min;
532 } // nsamples>0 check, some data found for this channel; not only trailer/header
533 } // end while over channel
535 }//end while over DDL's, of input stream
537 in->Reset(); // just in case the next customer forgets to check if the stream was reset..
539 // now check if it was a led event, only use high gain (that should be sufficient)
541 // by default all columns are accepted (init check to > 0)
542 int checkResultArray[AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALCols];
543 for (int ia=0; ia<(AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALCols); ia++) {
544 checkResultArray[ia] = 1;
547 if (fReqFractionAboveAmp) {
550 ok = CheckFractionAboveAmp(AmpValHighGain, checkResultArray);
552 if (!ok) return false; // skip event
555 fNAcceptedEvents++; // one more event accepted
557 if (fStartTime == 0) { // if start-timestamp wasn't set,we'll pick it up from the first event we encounter
558 fStartTime = Timestamp;
561 fHour = (Timestamp - fStartTime)/(double)fgkNumSecInHr;
562 if (fLatestHour < fHour) {
566 // it is a led event, now fill TTree
567 for(int i=0; i<fModules; i++){
568 for(int j=0; j<fColumns; j++){
569 if (checkResultArray[i*fColumns + j] > 0) { // column passed check
570 for(int k=0; k<fRows; k++){
572 TowerNum = GetTowerNum(i, j, k);
574 if(AmpValHighGain[TowerNum]) {
575 fAmp = AmpValHighGain[TowerNum];
576 fChannelNum = GetChannelNum(i,j,k,1);
577 fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValHighGain[TowerNum]);
578 fNHighGain[TowerNum]++;
580 if(AmpValLowGain[TowerNum]) {
581 fAmp = AmpValLowGain[TowerNum];
582 fChannelNum = GetChannelNum(i,j,k,0);
583 fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValLowGain[TowerNum]);
584 fNLowGain[TowerNum]++;
587 } // column passed check
590 // We also do the activity check for LEDRefs/Strips, but need to translate between column
591 // and strip indices for that; based on these relations:
592 // iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
593 // iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2;
595 // iColFirst = (iSM%2==0) ? iStrip*2 : (AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iStrip)*2;
598 for(int j=0; j<fLEDRefs; j++){
599 int iColFirst = (i%2==0) ? j*2 : (AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - j)*2; //CHECKME!!!
600 if ( (checkResultArray[i*fColumns + iColFirst]>0) || (checkResultArray[i*fColumns + iColFirst + 1]>0) ) { // at least one column in strip passed check
601 for (gain=0; gain<2; gain++) {
602 fRefNum = GetRefNum(i, j, gain);
603 if (LEDAmpVal[fRefNum]) {
604 fAmp = LEDAmpVal[fRefNum];
605 fTreeLEDAmpVsTime->Fill();//fRefNum,fHour,fAmp);
609 } // at least one column in strip passed check
617 //_____________________________________________________________________
618 Bool_t AliCaloCalibSignal::Save(TString fileName)
620 //Saves all the TTrees to the designated file
622 TFile destFile(fileName, "recreate");
624 if (destFile.IsZombie()) {
631 fTreeAmpVsTime->Write();
632 fTreeLEDAmpVsTime->Write();
634 Analyze(); // get the latest and greatest averages
635 fTreeAvgAmpVsTime->Write();
636 fTreeLEDAvgAmpVsTime->Write();
644 //_____________________________________________________________________
645 Bool_t AliCaloCalibSignal::Analyze()
647 // Fill the tree holding the average values
648 if (!fUseAverage) { return kFALSE; }
650 // Reset the average TTree if Analyze has already been called earlier,
651 // meaning that the TTree could have been partially filled
652 if (fTreeAvgAmpVsTime->GetEntries() > 0) {
653 fTreeAvgAmpVsTime->Reset();
656 //0: setup variables for the TProfile plots that we'll use to do the averages
660 if (fSecInAverage > 0) {
661 numProfBins = (int)( (fLatestHour*fgkNumSecInHr)/fSecInAverage + 1 ); // round-off
663 numProfBins += 2; // add extra buffer : first and last
664 double binSize = 1.0*fSecInAverage / fgkNumSecInHr;
666 timeMax = timeMin + numProfBins*binSize;
668 //1: set up TProfiles for the towers that had data
669 TProfile * profile[fgkMaxTowers*2]; // *2 is since we include both high and low gains
670 memset(profile, 0, sizeof(profile));
672 char name[200]; // for profile id and title
675 for (int i = 0; i<fModules; i++) {
676 for (int ic=0; ic<fColumns; ic++){
677 for (int ir=0; ir<fRows; ir++) {
679 TowerNum = GetTowerNum(i, ic, ir);
681 if (fNHighGain[TowerNum] > 0) {
682 fChannelNum = GetChannelNum(i, ic, ir, 1);
683 sprintf(name, "profileChan%d", fChannelNum);
684 profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
688 if (fNLowGain[TowerNum] > 0) {
689 fChannelNum = GetChannelNum(i, ic, ir, 0);
690 sprintf(name, "profileChan%d", fChannelNum);
691 profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
698 //2: fill profiles by looping over tree
699 // Set addresses for tree-readback also
700 fTreeAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
701 fTreeAmpVsTime->SetBranchAddress("fHour", &fHour);
702 fTreeAmpVsTime->SetBranchAddress("fAmp", &fAmp);
704 for (int ient=0; ient<fTreeAmpVsTime->GetEntries(); ient++) {
705 fTreeAmpVsTime->GetEntry(ient);
706 if (profile[fChannelNum]) {
707 // profile should always have been created above, for active channels
708 profile[fChannelNum]->Fill(fHour, fAmp);
712 // re-associating the branch addresses here seems to be needed for OK 'average' storage
713 fTreeAvgAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
714 fTreeAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
715 fTreeAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
716 fTreeAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);
718 //3: fill avg tree by looping over the profiles
719 for (fChannelNum = 0; fChannelNum<(fgkMaxTowers*2); fChannelNum++) {
720 if (profile[fChannelNum]) { // profile was created
721 if (profile[fChannelNum]->GetEntries() > 0) { // profile had some entries
722 for(int it=0; it<numProfBins; it++) {
723 if (profile[fChannelNum]->GetBinEntries(it+1) > 0) {
724 fAvgAmp = profile[fChannelNum]->GetBinContent(it+1);
725 fHour = profile[fChannelNum]->GetBinCenter(it+1);
726 fRMS = profile[fChannelNum]->GetBinError(it+1);
727 fTreeAvgAmpVsTime->Fill();
728 } // some entries for this bin
730 } // some entries for this profile
732 } // loop over all possible channels
735 // and finally, go through same exercise for LED also..
737 //1: set up TProfiles for the towers that had data
738 TProfile * profileLED[fgkMaxRefs*2]; // *2 is since we include both high and low gains
739 memset(profileLED, 0, sizeof(profileLED));
741 for (int i = 0; i<fModules; i++) {
742 for(int j=0; j<fLEDRefs; j++){
743 for (int gain=0; gain<2; gain++) {
744 fRefNum = GetRefNum(i, j, gain);
745 if (fNRef[fRefNum] > 0) {
746 sprintf(name, "profileLEDRef%d", fRefNum);
747 profileLED[fRefNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
753 //2: fill profiles by looping over tree
754 // Set addresses for tree-readback also
755 fTreeLEDAmpVsTime->SetBranchAddress("fRefNum", &fRefNum);
756 fTreeLEDAmpVsTime->SetBranchAddress("fHour", &fHour);
757 fTreeLEDAmpVsTime->SetBranchAddress("fAmp", &fAmp);
759 for (int ient=0; ient<fTreeLEDAmpVsTime->GetEntries(); ient++) {
760 fTreeLEDAmpVsTime->GetEntry(ient);
761 if (profileLED[fRefNum]) {
762 // profile should always have been created above, for active channels
763 profileLED[fRefNum]->Fill(fHour, fAmp);
767 // re-associating the branch addresses here seems to be needed for OK 'average' storage
768 fTreeLEDAvgAmpVsTime->SetBranchAddress("fRefNum", &fRefNum);
769 fTreeLEDAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
770 fTreeLEDAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
771 fTreeLEDAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);
773 //3: fill avg tree by looping over the profiles
774 for (fRefNum = 0; fRefNum<(fgkMaxRefs*2); fRefNum++) {
775 if (profileLED[fRefNum]) { // profile was created
776 if (profileLED[fRefNum]->GetEntries() > 0) { // profile had some entries
777 for(int it=0; it<numProfBins; it++) {
778 if (profileLED[fRefNum]->GetBinEntries(it+1) > 0) {
779 fAvgAmp = profileLED[fRefNum]->GetBinContent(it+1);
780 fHour = profileLED[fRefNum]->GetBinCenter(it+1);
781 fRMS = profileLED[fRefNum]->GetBinError(it+1);
782 fTreeLEDAvgAmpVsTime->Fill();
783 } // some entries for this bin
785 } // some entries for this profile
787 } // loop over all possible channels