Coverity corrections
[u/mrichter/AliRoot.git] / EMCAL / AliCaloCalibSignal.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
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 $ */
16
17 //________________________________________________________________________
18 //
19 // A help class for monitoring and calibration tools: MOOD, AMORE etc.,
20 // It can be created and used a la (ctor):
21 /*
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();
26 */
27 // fed an event:
28 //  fSignals->ProcessEvent(fCaloRawStream,fRawEventHeaderBase);
29 // get some info:
30 //  fSignals->GetXXX..()
31 // etc.
32 //________________________________________________________________________
33 #include <string>
34 #include <sstream>
35 #include <fstream>
36
37 #include "TProfile.h"
38 #include "TFile.h"
39
40 #include "AliRawReader.h"
41 #include "AliCaloRawStreamV3.h"
42
43 //The include file
44 #include "AliCaloCalibSignal.h"
45
46 ClassImp(AliCaloCalibSignal)
47
48 using namespace std;
49
50 // variables for TTree filling; not sure if they should be static or not
51 static int fChannelNum = 0; // for regular towers
52 static int fRefNum = 0; // for LED
53 static double fAmp = 0;
54 static double fAvgAmp = 0;
55 static double fRMS = 0;
56
57 // ctor; initialize everything in order to avoid compiler warnings
58 // put some reasonable defaults
59 AliCaloCalibSignal::AliCaloCalibSignal(kDetType detectorType) :  
60   TObject(),
61   fDetType(kNone),
62   fColumns(0),
63   fRows(0),
64   fLEDRefs(0),
65   fModules(0),
66   fCaloString(),
67   fMapping(NULL),
68   fRunNumber(-1),
69   fStartTime(0),
70   fAmpCut(40), // min. 40 ADC counts as default
71   fReqFractionAboveAmpCutVal(0.6), // 60% in a strip, per default
72   fReqFractionAboveAmp(kTRUE),
73   fAmpCutLEDRef(100), // min. 100 ADC counts as default
74   fReqLEDRefAboveAmpCutVal(kTRUE),
75   fHour(0),
76   fLatestHour(0),
77   fUseAverage(kTRUE),
78   fSecInAverage(1800), 
79   fNEvents(0),
80   fNAcceptedEvents(0),
81   fTreeAmpVsTime(NULL),
82   fTreeAvgAmpVsTime(NULL),
83   fTreeLEDAmpVsTime(NULL),
84   fTreeLEDAvgAmpVsTime(NULL)
85 {
86   //Default constructor. First we set the detector-type related constants.
87   if (detectorType == kPhos) {
88     fColumns = fgkPhosCols;
89     fRows = fgkPhosRows;
90     fLEDRefs = fgkPhosLEDRefs;
91     fModules = fgkPhosModules;
92     fCaloString = "PHOS";
93   } 
94   else {
95     //We'll just trust the enum to keep everything in line, so that if detectorType
96     //isn't kPhos then it is kEmCal. Note, however, that this is not necessarily the
97     //case, if someone intentionally gives another number
98     fColumns = AliEMCALGeoParams::fgkEMCALCols;
99     fRows = AliEMCALGeoParams::fgkEMCALRows;
100     fLEDRefs = AliEMCALGeoParams::fgkEMCALLEDRefs;
101     fModules = AliEMCALGeoParams::fgkEMCALModules;
102     fCaloString = "EMCAL";
103   }
104
105   fDetType = detectorType;
106
107   ResetInfo(); // trees and counters
108
109
110 // dtor
111 //_____________________________________________________________________
112 AliCaloCalibSignal::~AliCaloCalibSignal()
113 {
114   DeleteTrees();
115 }
116
117 //_____________________________________________________________________
118 void AliCaloCalibSignal::DeleteTrees()
119 {
120   // delete what was created in the ctor (TTrees)
121   if (fTreeAmpVsTime) delete fTreeAmpVsTime;
122   if (fTreeAvgAmpVsTime) delete fTreeAvgAmpVsTime;
123   if (fTreeLEDAmpVsTime) delete fTreeLEDAmpVsTime;
124   if (fTreeLEDAvgAmpVsTime) delete fTreeLEDAvgAmpVsTime;
125   // and reset pointers
126   fTreeAmpVsTime = NULL;
127   fTreeAvgAmpVsTime = NULL;
128   fTreeLEDAmpVsTime = NULL;
129   fTreeLEDAvgAmpVsTime = NULL;
130
131   return;
132 }
133
134 // copy ctor
135 //_____________________________________________________________________
136 AliCaloCalibSignal::AliCaloCalibSignal(const AliCaloCalibSignal &sig) :
137   TObject(sig),
138   fDetType(sig.GetDetectorType()),
139   fColumns(sig.GetColumns()),
140   fRows(sig.GetRows()),
141   fLEDRefs(sig.GetLEDRefs()),
142   fModules(sig.GetModules()),
143   fCaloString(sig.GetCaloString()),
144   fMapping(), //! note that we are not copying the map info
145   fRunNumber(sig.GetRunNumber()),
146   fStartTime(sig.GetStartTime()),
147   fAmpCut(sig.GetAmpCut()),
148   fReqFractionAboveAmpCutVal(sig.GetReqFractionAboveAmpCutVal()),
149   fReqFractionAboveAmp(sig.GetReqFractionAboveAmp()),
150   fAmpCutLEDRef(sig.GetAmpCutLEDRef()),
151   fReqLEDRefAboveAmpCutVal(sig.GetReqLEDRefAboveAmpCutVal()),
152   fHour(sig.GetHour()),
153   fLatestHour(sig.GetLatestHour()),
154   fUseAverage(sig.GetUseAverage()),
155   fSecInAverage(sig.GetSecInAverage()),
156   fNEvents(sig.GetNEvents()),
157   fNAcceptedEvents(sig.GetNAcceptedEvents()),
158   fTreeAmpVsTime(),
159   fTreeAvgAmpVsTime(),
160   fTreeLEDAmpVsTime(),
161   fTreeLEDAvgAmpVsTime()
162 {
163   // also the TTree contents
164   AddInfo(&sig);
165   for (Int_t i = 0; i<2*fgkMaxTowers; i++) {
166     if(i < fgkMaxTowers){
167       fNHighGain[i] = sig.fNHighGain[i];
168       fNLowGain[i]  = sig.fNLowGain[i]; 
169     }
170     fNRef[i] = sig.fNRef[i]; 
171   }
172   
173   
174   
175 }
176
177 // assignment operator; use copy ctor to make life easy..
178 //_____________________________________________________________________
179 AliCaloCalibSignal& AliCaloCalibSignal::operator = (const AliCaloCalibSignal &source)
180 {
181   // assignment operator; use copy ctor
182   if (&source == this) return *this;
183
184   new (this) AliCaloCalibSignal(source);
185   return *this;
186 }
187
188 //_____________________________________________________________________
189 void AliCaloCalibSignal::CreateTrees()
190 {
191   // initialize trees
192   // first, regular version
193   fTreeAmpVsTime = new TTree("fTreeAmpVsTime","Amplitude vs. Time Tree Variables");
194
195   fTreeAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I");
196   fTreeAmpVsTime->Branch("fHour", &fHour, "fHour/D");
197   fTreeAmpVsTime->Branch("fAmp", &fAmp, "fAmp/D");
198
199   // then, average version
200   fTreeAvgAmpVsTime = new TTree("fTreeAvgAmpVsTime","Average Amplitude vs. Time Tree Variables");
201
202   fTreeAvgAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I");
203   fTreeAvgAmpVsTime->Branch("fHour", &fHour, "fHour/D");
204   fTreeAvgAmpVsTime->Branch("fAvgAmp", &fAvgAmp, "fAvgAmp/D");
205   fTreeAvgAmpVsTime->Branch("fRMS", &fRMS, "fRMS/D");
206
207   // then same for LED..
208   fTreeLEDAmpVsTime = new TTree("fTreeLEDAmpVsTime","LED Amplitude vs. Time Tree Variables");
209   fTreeLEDAmpVsTime->Branch("fRefNum", &fRefNum, "fRefNum/I");
210   fTreeLEDAmpVsTime->Branch("fHour", &fHour, "fHour/D");
211   fTreeLEDAmpVsTime->Branch("fAmp", &fAmp, "fAmp/D");
212
213   fTreeLEDAvgAmpVsTime = new TTree("fTreeLEDAvgAmpVsTime","Average LED Amplitude vs. Time Tree Variables");
214   fTreeLEDAvgAmpVsTime->Branch("fRefNum", &fRefNum, "fRefNum/I");
215   fTreeLEDAvgAmpVsTime->Branch("fHour", &fHour, "fHour/D");
216   fTreeLEDAvgAmpVsTime->Branch("fAvgAmp", &fAvgAmp, "fAvgAmp/D");
217   fTreeLEDAvgAmpVsTime->Branch("fRMS", &fRMS, "fRMS/D");
218
219   return;
220 }
221
222 //_____________________________________________________________________
223 void AliCaloCalibSignal::ResetInfo()
224 { // reset trees and counters
225   Zero(); // set all counters to 0
226   DeleteTrees(); // delete previous stuff
227   CreateTrees(); // and create some new ones
228   return;
229 }
230
231 //_____________________________________________________________________
232 void AliCaloCalibSignal::Zero()
233 {
234   // set all counters to 0; not cuts etc. though
235   fHour = 0;
236   fLatestHour = 0;
237   fNEvents = 0;
238   fNAcceptedEvents = 0;
239
240   // Set the number of points for each tower: Amp vs. Time
241   memset(fNHighGain, 0, sizeof(fNHighGain));
242   memset(fNLowGain, 0, sizeof(fNLowGain));
243   // and LED reference
244   memset(fNRef, 0, sizeof(fNRef));
245
246   return;
247 }
248
249 //_____________________________________________________________________
250 Bool_t AliCaloCalibSignal::CheckFractionAboveAmp(const int *iAmpVal, 
251                                                  int resultArray[])
252 { // check fraction of towers, per column, that are above amplitude cut
253   Bool_t returnCode = false;
254     
255   int iTowerNum = 0;
256   double fraction = 0;
257   for (int i = 0; i<fModules; i++) {
258     for (int j = 0; j<fColumns; j++) {
259       int nAbove = 0;
260       for (int k = 0; k<fRows; k++) {
261         iTowerNum = GetTowerNum(i,j,k);
262         if (iAmpVal[iTowerNum] > fAmpCut) { 
263           nAbove++;
264         }
265       }
266       resultArray[i*fColumns +j] = 0; // init. to denied
267       if (nAbove > 0) {
268         fraction = (1.0*nAbove) / fRows;
269         /*
270         printf("DS mod %d col %d nAbove %d fraction %3.2f\n",
271                i, j, nAbove, fraction);
272         */
273         if (fraction > fReqFractionAboveAmpCutVal) {
274           resultArray[i*fColumns + j] = nAbove;
275           returnCode = true;
276         }  
277       }
278     }
279   } // modules loop
280   
281   return returnCode;
282 }
283
284
285 //_____________________________________________________________________
286 Bool_t AliCaloCalibSignal::CheckLEDRefAboveAmp(const int *iAmpVal, 
287                                                int resultArray[])
288 { // check which LEDRef/Mon strips are above amplitude cut
289   Bool_t returnCode = false;
290     
291   int iRefNum = 0;
292   int gain = 1; // look at high gain; this should be rather saturated usually..
293   for (int i = 0; i<fModules; i++) {
294     for (int j = 0; j<fLEDRefs; j++) {
295       iRefNum = GetRefNum(i, j, gain);
296       if (iAmpVal[iRefNum] > fAmpCutLEDRef) { 
297         resultArray[i*fLEDRefs +j] = 1; // enough signal
298         returnCode = true;
299       }
300       else {
301         resultArray[i*fLEDRefs +j] = 0; // not enough signal
302       }
303       
304       /*
305       printf("DS mod %d LEDRef %d ampVal %d\n",
306              i, j, iAmpVal[iRefNum]);
307       */
308     } // LEDRefs
309   } // modules loop
310   
311   return returnCode;
312 }
313
314 // Parameter/cut handling
315 //_____________________________________________________________________
316 void AliCaloCalibSignal::SetParametersFromFile(const char *parameterFile)
317 { // set parameters from file
318   static const string delimitor("::");
319         
320   // open, check input file
321   ifstream in( parameterFile );
322   if( !in ) {
323     printf("in AliCaloCalibSignal::SetParametersFromFile - Using default/run_time parameters.\n");
324     return;
325   } 
326
327   // Note: this method is a bit more complicated than it really has to be
328   // - allowing for multiple entries per line, arbitrary order of the
329   // different variables etc. But I wanted to try and do this in as
330   // correct a C++ way as I could (as an exercise).
331
332   // read in
333   char readline[1024];
334   while ((in.rdstate() & ios::failbit) == 0 ) {
335     
336     // Read into the raw char array and then construct a string
337     // to do the searching
338     in.getline(readline, 1024);
339     istringstream s(readline);          
340                 
341     while ( ( s.rdstate() & ios::failbit ) == 0 ) {
342                         
343       string keyValue; 
344       s >> keyValue;
345       
346       // check stream status
347       if( ( s.rdstate() & ios::failbit ) == ios::failbit ) break;
348                         
349       // skip rest of line if comments found
350       if( keyValue.substr( 0, 2 ) == "//" ) break;
351                         
352       // look for "::" in keyValue pair
353       size_t position = keyValue.find( delimitor );
354       if( position == string::npos ) {
355         printf("wrong format for key::value pair: %s\n", keyValue.c_str());
356       }
357                                 
358       // split keyValue pair
359       string key( keyValue.substr( 0, position ) );
360       string value( keyValue.substr( position+delimitor.size(), 
361                                       keyValue.size()-delimitor.size() ) );
362                         
363       // check value does not contain a new delimitor
364       if( value.find( delimitor ) != string::npos ) {
365         printf("wrong format for key::value pair: %s\n", keyValue.c_str());
366       }
367       
368       // debug: check key value pair
369       // printf("AliCaloCalibSignal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
370
371       // if the key matches with something we expect, we assign the new value
372       if ( (key == "fAmpCut") || (key == "fReqFractionAboveAmpCutVal") ||
373            (key == "fAmpCutLEDRef") || (key == "fSecInAverage") ) {
374         istringstream iss(value);
375         printf("AliCaloCalibSignal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
376
377         if (key == "fAmpCut") { 
378           iss >> fAmpCut; 
379         }
380         else if (key == "fReqFractionAboveAmpCutVal") { 
381           iss >> fReqFractionAboveAmpCutVal; 
382         }
383         else if (key == "fAmpCutLEDRef") { 
384           iss >> fAmpCutLEDRef; 
385         }
386         else if (key == "fSecInAverage") { 
387           iss >> fSecInAverage; 
388         }
389       } // some match found/expected
390
391     }           
392   }
393
394   in.close();
395   return;       
396 }
397
398 //_____________________________________________________________________
399 void AliCaloCalibSignal::WriteParametersToFile(const char *parameterFile)
400 { // write parameters to file
401   static const string delimitor("::");
402   ofstream out( parameterFile );
403   out << "// " << parameterFile << endl;
404   out << "fAmpCut" << "::" << fAmpCut << endl;
405   out << "fReqFractionAboveAmpCutVal" << "::" << fReqFractionAboveAmpCutVal << endl;
406   out << "fAmpCutLEDRef" << "::" << fAmpCutLEDRef << endl;
407   out << "fSecInAverage" << "::" << fSecInAverage << endl;
408
409   out.close();
410   return;
411 }
412
413 //_____________________________________________________________________
414 Bool_t AliCaloCalibSignal::AddInfo(const AliCaloCalibSignal *sig)
415 {
416   // 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.
417
418   // add info from sig's TTrees to ours..
419   TTree *sigAmp = sig->GetTreeAmpVsTime();
420   TTree *sigAvgAmp = sig->GetTreeAvgAmpVsTime();
421
422   // we could try some merging via TList or what also as a more elegant approach
423   // but I wanted with the stupid/simple and hopefully safe approach of looping
424   // over what we want to add..
425
426   // associate variables for sigAmp and sigAvgAmp:
427   sigAmp->SetBranchAddress("fChannelNum",&fChannelNum);
428   sigAmp->SetBranchAddress("fHour",&fHour);
429   sigAmp->SetBranchAddress("fAmp",&fAmp);
430
431   // loop over the trees.. note that since we use the same variables we should not need
432   // to do any assignments between the getting and filling
433   for (int i=0; i<sigAmp->GetEntries(); i++) {
434     sigAmp->GetEntry(i);
435     fTreeAmpVsTime->Fill();
436   }
437
438   sigAvgAmp->SetBranchAddress("fChannelNum",&fChannelNum);
439   sigAvgAmp->SetBranchAddress("fHour",&fHour);
440   sigAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
441   sigAvgAmp->SetBranchAddress("fRMS",&fRMS);
442
443   for (int i=0; i<sigAvgAmp->GetEntries(); i++) {
444     sigAvgAmp->GetEntry(i);
445     fTreeAvgAmpVsTime->Fill();
446   }
447
448   // also LED.. 
449   TTree *sigLEDAmp = sig->GetTreeLEDAmpVsTime();
450   TTree *sigLEDAvgAmp = sig->GetTreeLEDAvgAmpVsTime();
451
452   // associate variables for sigAmp and sigAvgAmp:
453   sigLEDAmp->SetBranchAddress("fRefNum",&fRefNum);
454   sigLEDAmp->SetBranchAddress("fHour",&fHour);
455   sigLEDAmp->SetBranchAddress("fAmp",&fAmp);
456
457   // loop over the trees.. note that since we use the same variables we should not need
458   // to do any assignments between the getting and filling
459   for (int i=0; i<sigLEDAmp->GetEntries(); i++) {
460     sigLEDAmp->GetEntry(i);
461     fTreeLEDAmpVsTime->Fill();
462   }
463
464   sigLEDAvgAmp->SetBranchAddress("fRefNum",&fRefNum);
465   sigLEDAvgAmp->SetBranchAddress("fHour",&fHour);
466   sigLEDAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
467   sigLEDAvgAmp->SetBranchAddress("fRMS",&fRMS);
468
469   for (int i=0; i<sigLEDAvgAmp->GetEntries(); i++) {
470     sigLEDAvgAmp->GetEntry(i);
471     fTreeLEDAvgAmpVsTime->Fill();
472   }
473
474
475   return kTRUE;//We hopefully succesfully added info from the supplied object
476 }
477
478 //_____________________________________________________________________
479 Bool_t AliCaloCalibSignal::ProcessEvent(AliRawReader *rawReader)
480 {
481   // if fMapping is NULL the rawstream will crate its own mapping
482   AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);  
483   if (fDetType == kEmCal) {
484     rawReader->Select("EMCAL", 0, AliEMCALGeoParams::fgkLastAltroDDL) ; //select EMCAL DDL range 
485   }
486   return ProcessEvent( &rawStream, rawReader->GetTimestamp() );
487 }
488
489 //_____________________________________________________________________
490 Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStreamV3 *in, UInt_t Timestamp)
491
492   // Method to process=analyze one event in the data stream
493   if (!in) return kFALSE; //Return right away if there's a null pointer
494   
495   fNEvents++; // one more event
496
497   // use maximum numbers to set array sizes
498   int iAmpValHighGain[fgkMaxTowers];
499   int iAmpValLowGain[fgkMaxTowers];
500   memset(iAmpValHighGain, 0, sizeof(iAmpValHighGain));
501   memset(iAmpValLowGain, 0, sizeof(iAmpValLowGain));
502
503   // also for LED reference
504   int iLEDAmpVal[fgkMaxRefs * 2]; // factor 2 is for the two gain values
505   memset(iLEDAmpVal, 0, sizeof(iLEDAmpVal));
506
507   int sample = 0; // temporary value
508   int gain = 0; // high or low gain
509   
510   // Number of Low and High gain, and LED Ref, channels for this event:
511   int nLowChan = 0; 
512   int nHighChan = 0; 
513   int nLEDRefChan = 0;
514
515   int iTowerNum = 0; // array index for regular towers
516   int iRefNum = 0; // array index for LED references
517
518   // loop first to get the fraction of channels with amplitudes above cut
519
520   while (in->NextDDL()) {
521     while (in->NextChannel()) {
522
523       // counters
524       int max = AliEMCALGeoParams::fgkSampleMin, min = AliEMCALGeoParams::fgkSampleMax; // min and max sample values
525       int nsamples = 0;
526
527       while (in->NextBunch()) {
528         const UShort_t *sig = in->GetSignals();
529         nsamples += in->GetBunchLength();
530         for (Int_t i = 0; i < in->GetBunchLength(); i++) {
531           sample = sig[i];
532
533           // check if it's a min or max value
534           if (sample < min) min = sample;
535           if (sample > max) max = sample;
536
537         } // loop over samples in bunch
538       } // loop over bunches
539
540       if (nsamples > 0) { // this check is needed for when we have zero-supp. on, but not sparse readout
541
542       gain = -1; // init to not valid value
543       //If we're here then we're done with this tower
544       if ( in->IsLowGain() ) {
545         gain = 0;
546       }
547       else if ( in->IsHighGain() ) {
548         gain = 1;
549       }
550       else if ( in->IsLEDMonData() ) {
551         gain = in->GetRow(); // gain coded in (in RCU/Altro mapping) as Row info for LED refs..
552       }
553
554       // it should be enough to check the SuperModule info for each DDL really, but let's keep it here for now
555       int arrayPos = in->GetModule(); //The modules are numbered starting from 0
556       //Debug
557       if (arrayPos < 0 || arrayPos >= fModules) {
558         printf("AliCaloCalibSignal::ProcessEvent = Oh no: arrayPos = %i.\n", arrayPos); 
559         return kFALSE;
560       }
561
562       if ( in->IsHighGain() || in->IsLowGain() ) { // regular tower
563         // get tower number for AmpVal array
564         iTowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow()); 
565
566         if (gain == 0) {
567           // fill amplitude into the array         
568           iAmpValLowGain[iTowerNum] = max - min;
569           nLowChan++;
570         } 
571         else if (gain==1) {//fill the high gain ones
572           // fill amplitude into the array
573           iAmpValHighGain[iTowerNum] = max - min;
574           nHighChan++;
575         }//end if gain
576       } // regular tower
577       else if ( in->IsLEDMonData() ) { // LED ref.; 
578         // strip # is coded is 'column' in the channel maps 
579         iRefNum = GetRefNum(arrayPos, in->GetColumn(), gain); 
580         iLEDAmpVal[iRefNum] = max - min;
581         nLEDRefChan++;
582       } // end of LED ref
583
584       } // nsamples>0 check, some data found for this channel; not only trailer/header      
585     } // end while over channel 
586    
587   }//end while over DDL's, of input stream
588   
589   in->Reset(); // just in case the next customer forgets to check if the stream was reset..
590
591   // now check if it was an LED event, using the LED Reference info per strip
592
593   // by default all columns are accepted (init check to > 0)
594   int checkResultArray[AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALCols];
595   for (int ia=0; ia<(AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALCols); ia++) { 
596     checkResultArray[ia] = 1; 
597   }
598   if (fReqFractionAboveAmp) {
599     bool ok = false;
600     if (nHighChan > 0) { 
601       ok = CheckFractionAboveAmp(iAmpValHighGain, checkResultArray); 
602     }
603     if (!ok) return false; // skip event
604   }
605
606   // by default all columns are accepted (init check to > 0)
607   int checkResultArrayLEDRef[AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALLEDRefs];
608   for (int ia=0; ia<(AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALLEDRefs); ia++) { 
609     checkResultArrayLEDRef[ia] = 1; 
610   }
611   if (fReqLEDRefAboveAmpCutVal) {
612     bool ok = false;
613     if (nLEDRefChan > 0) { 
614       ok = CheckLEDRefAboveAmp(iLEDAmpVal, checkResultArrayLEDRef); 
615     }
616     if (!ok) return false; // skip event
617   }
618
619   fNAcceptedEvents++; // one more event accepted
620
621   if (fStartTime == 0) { // if start-timestamp wasn't set,we'll pick it up from the first event we encounter
622     fStartTime = Timestamp;
623   }
624
625   fHour = (Timestamp - fStartTime)/(double)fgkNumSecInHr;
626   if (fLatestHour < fHour) {
627     fLatestHour = fHour; 
628   }
629   
630   // it is a led event, now fill TTree
631   // We also do the activity check for LEDRefs/Strips, but need to translate between column
632   // and strip indices for that; based on these relations: 
633   // iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
634   // iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2;
635   // which leads to
636   // iColFirst = (iSM%2==0) ? iStrip*2 : (AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iStrip)*2;
637
638   for(int i=0; i<fModules; i++){
639     for(int j=0; j<fColumns; j++) {
640       int iStrip = (i%2==0) ? j/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - j/2;
641       if (checkResultArray[i*fColumns + j]>0  && checkResultArrayLEDRef[i*fLEDRefs + iStrip]>0) { // column passed check 
642       for(int k=0; k<fRows; k++){
643         
644         iTowerNum = GetTowerNum(i, j, k); 
645
646         if(iAmpValHighGain[iTowerNum]) {
647           fAmp = iAmpValHighGain[iTowerNum];
648           fChannelNum = GetChannelNum(i,j,k,1);
649           fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValHighGain[iTowerNum]);
650           fNHighGain[iTowerNum]++;
651         }
652         if(iAmpValLowGain[iTowerNum]) {
653           fAmp = iAmpValLowGain[iTowerNum];
654           fChannelNum = GetChannelNum(i,j,k,0);
655           fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValLowGain[iTowerNum]);
656           fNLowGain[iTowerNum]++;
657         }
658       } // rows
659       } // column passed check, and LED Ref for strip passed check (if any)
660     } // columns
661
662     // also LED refs
663     for(int j=0; j<fLEDRefs; j++){
664       int iColFirst = (i%2==0) ? j*2 : (AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - j)*2; //CHECKME!!!
665       if ( ((checkResultArray[i*fColumns + iColFirst]>0) || (checkResultArray[i*fColumns + iColFirst + 1]>0)) && // at least one column in strip passed check 
666            (checkResultArrayLEDRef[i*fLEDRefs + j]>0) ) { // and LED Ref passed checks
667         for (gain=0; gain<2; gain++) {
668           fRefNum = GetRefNum(i, j, gain); 
669           if (iLEDAmpVal[fRefNum]) {
670             fAmp = iLEDAmpVal[fRefNum];
671             fTreeLEDAmpVsTime->Fill();//fRefNum,fHour,fAmp);
672             fNRef[fRefNum]++;
673           }
674         } // gain
675       } // at least one column in strip passed check, and LED Ref passed check (if any) 
676     }
677
678   } // modules
679   
680   return kTRUE;
681 }
682
683 //_____________________________________________________________________
684 Bool_t AliCaloCalibSignal::Save(TString fileName)
685 {
686   //Saves all the TTrees to the designated file
687   
688   TFile destFile(fileName, "recreate");
689   
690   if (destFile.IsZombie()) {
691     return kFALSE;
692   }
693   
694   destFile.cd();
695
696   // save the trees
697   fTreeAmpVsTime->Write();
698   fTreeLEDAmpVsTime->Write();
699   if (fUseAverage) { 
700     Analyze(); // get the latest and greatest averages
701     fTreeAvgAmpVsTime->Write();
702     fTreeLEDAvgAmpVsTime->Write();
703   }
704
705   destFile.Close();
706   
707   return kTRUE;
708 }
709
710 //_____________________________________________________________________
711 Bool_t AliCaloCalibSignal::Analyze()
712 {
713   // Fill the tree holding the average values
714   if (!fUseAverage) { return kFALSE; }
715
716   // Reset the average TTree if Analyze has already been called earlier,
717   // meaning that the TTree could have been partially filled
718   if (fTreeAvgAmpVsTime->GetEntries() > 0) {
719     fTreeAvgAmpVsTime->Reset();
720   }
721
722   //0: setup variables for the TProfile plots that we'll use to do the averages
723   int numProfBins = 0;
724   double timeMin = 0;
725   double timeMax = 0;
726   if (fSecInAverage > 0) {
727     numProfBins = (int)( (fLatestHour*fgkNumSecInHr)/fSecInAverage + 1 ); // round-off
728   }
729   numProfBins += 2; // add extra buffer : first and last
730   double binSize = 1.0*fSecInAverage / fgkNumSecInHr;
731   timeMin = - binSize;
732   timeMax = timeMin + numProfBins*binSize;
733
734   //1: set up TProfiles for the towers that had data
735   TProfile * profile[fgkMaxTowers*2]; // *2 is since we include both high and low gains
736   memset(profile, 0, sizeof(profile));
737   const Int_t buffersize = 200;
738   char name[buffersize]; // for profile id and title
739   int iTowerNum = 0;
740
741   for (int i = 0; i<fModules; i++) {
742     for (int ic=0; ic<fColumns; ic++){
743       for (int ir=0; ir<fRows; ir++) {
744
745         iTowerNum = GetTowerNum(i, ic, ir);
746         // high gain
747         if (fNHighGain[iTowerNum] > 0) {
748           fChannelNum = GetChannelNum(i, ic, ir, 1); 
749           snprintf(name,buffersize,"profileChan%d", fChannelNum);
750           profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
751         }
752
753         // same for low gain
754         if (fNLowGain[iTowerNum] > 0) {
755           fChannelNum = GetChannelNum(i, ic, ir, 0); 
756           snprintf(name,buffersize,"profileChan%d", fChannelNum);
757           profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
758         }
759
760       } // rows
761     } // columns
762   } // modules
763
764   //2: fill profiles by looping over tree
765   // Set addresses for tree-readback also
766   fTreeAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
767   fTreeAmpVsTime->SetBranchAddress("fHour", &fHour);
768   fTreeAmpVsTime->SetBranchAddress("fAmp", &fAmp);
769
770   for (int ient=0; ient<fTreeAmpVsTime->GetEntries(); ient++) {
771     fTreeAmpVsTime->GetEntry(ient);
772     if (profile[fChannelNum]) { 
773       // profile should always have been created above, for active channels
774       profile[fChannelNum]->Fill(fHour, fAmp);
775     }
776   }
777
778   // re-associating the branch addresses here seems to be needed for OK 'average' storage           
779   fTreeAvgAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
780   fTreeAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
781   fTreeAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
782   fTreeAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);
783
784   //3: fill avg tree by looping over the profiles
785   for (fChannelNum = 0; fChannelNum<(fgkMaxTowers*2); fChannelNum++) {
786     if (profile[fChannelNum]) { // profile was created
787       if (profile[fChannelNum]->GetEntries() > 0) { // profile had some entries
788         for(int it=0; it<numProfBins; it++) {
789           if (profile[fChannelNum]->GetBinEntries(it+1) > 0) {
790             fAvgAmp = profile[fChannelNum]->GetBinContent(it+1);
791             fHour = profile[fChannelNum]->GetBinCenter(it+1);
792             fRMS = profile[fChannelNum]->GetBinError(it+1);
793             fTreeAvgAmpVsTime->Fill();
794           } // some entries for this bin
795         } // loop over bins
796       } // some entries for this profile
797     } // profile exists  
798   } // loop over all possible channels
799
800
801   // and finally, go through same exercise for LED also.. 
802
803   //1: set up TProfiles for the towers that had data
804   TProfile * profileLED[fgkMaxRefs*2]; // *2 is since we include both high and low gains
805   memset(profileLED, 0, sizeof(profileLED));
806
807   for (int i = 0; i<fModules; i++) {
808     for(int j=0; j<fLEDRefs; j++){
809       for (int gain=0; gain<2; gain++) {
810         fRefNum = GetRefNum(i, j, gain);
811         if (fNRef[fRefNum] > 0) { 
812           sprintf(name, "profileLEDRef%d", fRefNum);
813           profileLED[fRefNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
814         } 
815       }// gain
816     } 
817   } // modules
818
819   //2: fill profiles by looping over tree
820   // Set addresses for tree-readback also
821   fTreeLEDAmpVsTime->SetBranchAddress("fRefNum", &fRefNum);
822   fTreeLEDAmpVsTime->SetBranchAddress("fHour", &fHour);
823   fTreeLEDAmpVsTime->SetBranchAddress("fAmp", &fAmp);
824
825   for (int ient=0; ient<fTreeLEDAmpVsTime->GetEntries(); ient++) {
826     fTreeLEDAmpVsTime->GetEntry(ient);
827     if (profileLED[fRefNum]) { 
828       // profile should always have been created above, for active channels
829       profileLED[fRefNum]->Fill(fHour, fAmp);
830     }
831   }
832
833   // re-associating the branch addresses here seems to be needed for OK 'average' storage           
834   fTreeLEDAvgAmpVsTime->SetBranchAddress("fRefNum", &fRefNum);
835   fTreeLEDAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
836   fTreeLEDAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
837   fTreeLEDAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);
838
839   //3: fill avg tree by looping over the profiles
840   for (fRefNum = 0; fRefNum<(fgkMaxRefs*2); fRefNum++) {
841     if (profileLED[fRefNum]) { // profile was created
842       if (profileLED[fRefNum]->GetEntries() > 0) { // profile had some entries
843         for(int it=0; it<numProfBins; it++) {
844           if (profileLED[fRefNum]->GetBinEntries(it+1) > 0) {
845             fAvgAmp = profileLED[fRefNum]->GetBinContent(it+1);
846             fHour = profileLED[fRefNum]->GetBinCenter(it+1);
847             fRMS = profileLED[fRefNum]->GetBinError(it+1);
848             fTreeLEDAvgAmpVsTime->Fill();
849           } // some entries for this bin
850         } // loop over bins
851       } // some entries for this profile
852     } // profile exists  
853   } // loop over all possible channels
854
855   // OK, we're done..
856
857   return kTRUE;
858 }
859
860 //_____________________________________________________________________
861 Bool_t AliCaloCalibSignal::DecodeChannelNum(const int chanId, 
862                                             int *imod, int *icol, int *irow, int *igain) const  
863 { // return the module, column, row, and gain for a given channel number
864   *igain = chanId/(fModules*fColumns*fRows);
865   *imod = (chanId/(fColumns*fRows)) % fModules;
866   *icol = (chanId/fRows) % fColumns;
867   *irow = chanId % fRows;
868   return kTRUE;
869
870
871 //_____________________________________________________________________
872 Bool_t AliCaloCalibSignal::DecodeRefNum(const int refId, 
873                                         int *imod, int *istripMod, int *igain) const 
874 { // return the module, stripModule, and gain for a given reference number
875   *igain = refId/(fModules*fLEDRefs);
876   *imod = (refId/(fLEDRefs)) % fModules;
877   *istripMod = refId % fLEDRefs;
878   return kTRUE;
879