bug fix in the vertex selection
[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; // for regular towers
52 static int fRefNum; // for LED
53 static double fAmp;
54 static double fAvgAmp;
55 static double fRMS;
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   fHour(0),
74   fLatestHour(0),
75   fUseAverage(kTRUE),
76   fSecInAverage(1800), 
77   fNEvents(0),
78   fNAcceptedEvents(0),
79   fTreeAmpVsTime(NULL),
80   fTreeAvgAmpVsTime(NULL),
81   fTreeLEDAmpVsTime(NULL),
82   fTreeLEDAvgAmpVsTime(NULL)
83 {
84   //Default constructor. First we set the detector-type related constants.
85   if (detectorType == kPhos) {
86     fColumns = fgkPhosCols;
87     fRows = fgkPhosRows;
88     fLEDRefs = fgkPhosLEDRefs;
89     fModules = fgkPhosModules;
90     fCaloString = "PHOS";
91   } 
92   else {
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";
101   }
102
103   fDetType = detectorType;
104
105   ResetInfo(); // trees and counters
106
107
108 // dtor
109 //_____________________________________________________________________
110 AliCaloCalibSignal::~AliCaloCalibSignal()
111 {
112   DeleteTrees();
113 }
114
115 //_____________________________________________________________________
116 void AliCaloCalibSignal::DeleteTrees()
117 {
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;
128
129   return;
130 }
131
132 // copy ctor
133 //_____________________________________________________________________
134 AliCaloCalibSignal::AliCaloCalibSignal(const AliCaloCalibSignal &sig) :
135   TObject(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)
158 {
159   // also the TTree contents
160   AddInfo(&sig);
161 }
162
163 // assignment operator; use copy ctor to make life easy..
164 //_____________________________________________________________________
165 AliCaloCalibSignal& AliCaloCalibSignal::operator = (const AliCaloCalibSignal &source)
166 {
167   // assignment operator; use copy ctor
168   if (&source == this) return *this;
169
170   new (this) AliCaloCalibSignal(source);
171   return *this;
172 }
173
174 //_____________________________________________________________________
175 void AliCaloCalibSignal::CreateTrees()
176 {
177   // initialize trees
178   // first, regular version
179   fTreeAmpVsTime = new TTree("fTreeAmpVsTime","Amplitude vs. Time Tree Variables");
180
181   fTreeAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I");
182   fTreeAmpVsTime->Branch("fHour", &fHour, "fHour/D");
183   fTreeAmpVsTime->Branch("fAmp", &fAmp, "fAmp/D");
184
185   // then, average version
186   fTreeAvgAmpVsTime = new TTree("fTreeAvgAmpVsTime","Average Amplitude vs. Time Tree Variables");
187
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");
192
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");
198
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");
204
205   return;
206 }
207
208 //_____________________________________________________________________
209 void AliCaloCalibSignal::ResetInfo()
210 { // reset trees and counters
211   Zero(); // set all counters to 0
212   DeleteTrees(); // delete previous stuff
213   CreateTrees(); // and create some new ones
214   return;
215 }
216
217 //_____________________________________________________________________
218 void AliCaloCalibSignal::Zero()
219 {
220   // set all counters to 0; not cuts etc. though
221   fHour = 0;
222   fLatestHour = 0;
223   fNEvents = 0;
224   fNAcceptedEvents = 0;
225
226   // Set the number of points for each tower: Amp vs. Time
227   memset(fNHighGain, 0, sizeof(fNHighGain));
228   memset(fNLowGain, 0, sizeof(fNLowGain));
229   // and LED reference
230   memset(fNRef, 0, sizeof(fNRef));
231
232   return;
233 }
234
235 //_____________________________________________________________________
236 Bool_t AliCaloCalibSignal::CheckFractionAboveAmp(const int *iAmpVal, 
237                                                  int resultArray[])
238 { // check fraction of towers, per column, that are above amplitude cut
239   Bool_t returnCode = false;
240     
241   int iTowerNum = 0;
242   double fraction = 0;
243   for (int i = 0; i<fModules; i++) {
244     for (int j = 0; j<fColumns; j++) {
245       int nAbove = 0;
246       for (int k = 0; k<fRows; k++) {
247         iTowerNum = GetTowerNum(i,j,k);
248         if (iAmpVal[iTowerNum] > fAmpCut) { 
249           nAbove++;
250         }
251       }
252       resultArray[i*fColumns +j] = 0; // init. to denied
253       if (nAbove > 0) {
254         fraction = (1.0*nAbove) / fRows;
255         /*
256         printf("DS mod %d col %d nAbove %d fraction %3.2f\n",
257                i, j, nAbove, fraction);
258         */
259         if (fraction > fReqFractionAboveAmpCutVal) {
260           resultArray[i*fColumns + j] = nAbove;
261           returnCode = true;
262         }  
263       }
264     }
265   } // modules loop
266   
267   return returnCode;
268 }
269
270 // Parameter/cut handling
271 //_____________________________________________________________________
272 void AliCaloCalibSignal::SetParametersFromFile(const char *parameterFile)
273 { // set parameters from file
274   static const string delimitor("::");
275         
276   // open, check input file
277   ifstream in( parameterFile );
278   if( !in ) {
279     printf("in AliCaloCalibSignal::SetParametersFromFile - Using default/run_time parameters.\n");
280     return;
281   } 
282
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).
287
288   // read in
289   char readline[1024];
290   while ((in.rdstate() & ios::failbit) == 0 ) {
291     
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);          
296                 
297     while ( ( s.rdstate() & ios::failbit ) == 0 ) {
298                         
299       string keyValue; 
300       s >> keyValue;
301       
302       // check stream status
303       if( s.rdstate() & ios::failbit ) break;
304                         
305       // skip rest of line if comments found
306       if( keyValue.substr( 0, 2 ) == "//" ) break;
307                         
308       // look for "::" in keyValue pair
309       size_t position = keyValue.find( delimitor );
310       if( position == string::npos ) {
311         printf("wrong format for key::value pair: %s\n", keyValue.c_str());
312       }
313                                 
314       // split keyValue pair
315       string key( keyValue.substr( 0, position ) );
316       string value( keyValue.substr( position+delimitor.size(), 
317                                       keyValue.size()-delimitor.size() ) );
318                         
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", keyValue.c_str());
322       }
323       
324       // debug: check key value pair
325       // printf("AliCaloCalibSignal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
326
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());
331
332         if (key == "fAmpCut") { 
333           iss >> fAmpCut; 
334         }
335         else if (key == "fReqFractionAboveAmpCutVal") { 
336           iss >> fReqFractionAboveAmpCutVal; 
337         }
338         else if (key == "fSecInAverage") { 
339           iss >> fSecInAverage; 
340         }
341       } // some match found/expected
342
343     }           
344   }
345
346   in.close();
347   return;
348         
349 }
350
351 //_____________________________________________________________________
352 void AliCaloCalibSignal::WriteParametersToFile(const char *parameterFile)
353 { // write parameters to file
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;
360
361   out.close();
362   return;
363 }
364
365 //_____________________________________________________________________
366 Bool_t AliCaloCalibSignal::AddInfo(const AliCaloCalibSignal *sig)
367 {
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.
369
370   // add info from sig's TTrees to ours..
371   TTree *sigAmp = sig->GetTreeAmpVsTime();
372   TTree *sigAvgAmp = sig->GetTreeAvgAmpVsTime();
373
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..
377
378   // associate variables for sigAmp and sigAvgAmp:
379   sigAmp->SetBranchAddress("fChannelNum",&fChannelNum);
380   sigAmp->SetBranchAddress("fHour",&fHour);
381   sigAmp->SetBranchAddress("fAmp",&fAmp);
382
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++) {
386     sigAmp->GetEntry(i);
387     fTreeAmpVsTime->Fill();
388   }
389
390   sigAvgAmp->SetBranchAddress("fChannelNum",&fChannelNum);
391   sigAvgAmp->SetBranchAddress("fHour",&fHour);
392   sigAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
393   sigAvgAmp->SetBranchAddress("fRMS",&fRMS);
394
395   for (int i=0; i<sigAvgAmp->GetEntries(); i++) {
396     sigAvgAmp->GetEntry(i);
397     fTreeAvgAmpVsTime->Fill();
398   }
399
400   // also LED.. 
401   TTree *sigLEDAmp = sig->GetTreeLEDAmpVsTime();
402   TTree *sigLEDAvgAmp = sig->GetTreeLEDAvgAmpVsTime();
403
404   // associate variables for sigAmp and sigAvgAmp:
405   sigLEDAmp->SetBranchAddress("fRefNum",&fRefNum);
406   sigLEDAmp->SetBranchAddress("fHour",&fHour);
407   sigLEDAmp->SetBranchAddress("fAmp",&fAmp);
408
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();
414   }
415
416   sigLEDAvgAmp->SetBranchAddress("fRefNum",&fRefNum);
417   sigLEDAvgAmp->SetBranchAddress("fHour",&fHour);
418   sigLEDAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
419   sigLEDAvgAmp->SetBranchAddress("fRMS",&fRMS);
420
421   for (int i=0; i<sigLEDAvgAmp->GetEntries(); i++) {
422     sigLEDAvgAmp->GetEntry(i);
423     fTreeLEDAvgAmpVsTime->Fill();
424   }
425
426
427   return kTRUE;//We hopefully succesfully added info from the supplied object
428 }
429
430 //_____________________________________________________________________
431 Bool_t AliCaloCalibSignal::ProcessEvent(AliRawReader *rawReader)
432 {
433   // if fMapping is NULL the rawstream will crate its own mapping
434   AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
435
436   return ProcessEvent( &rawStream, rawReader->GetTimestamp() );
437 }
438
439 //_____________________________________________________________________
440 Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStreamV3 *in, UInt_t Timestamp)
441
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
444   
445   fNEvents++; // one more event
446
447   // use maximum numbers to set array sizes
448   int iAmpValHighGain[fgkMaxTowers];
449   int iAmpValLowGain[fgkMaxTowers];
450   memset(iAmpValHighGain, 0, sizeof(iAmpValHighGain));
451   memset(iAmpValLowGain, 0, sizeof(iAmpValLowGain));
452
453   // also for LED reference
454   int iLEDAmpVal[fgkMaxRefs * 2]; // factor 2 is for the two gain values
455   memset(iLEDAmpVal, 0, sizeof(iLEDAmpVal));
456
457   int sample; // temporary value
458   int gain = 0; // high or low gain
459   
460   // Number of Low and High gain channels for this event:
461   int nLowChan = 0; 
462   int nHighChan = 0; 
463
464   int iTowerNum = 0; // array index for regular towers
465   int iRefNum = 0; // array index for LED references
466
467   // loop first to get the fraction of channels with amplitudes above cut
468
469   while (in->NextDDL()) {
470     while (in->NextChannel()) {
471
472       // counters
473       int max = AliEMCALGeoParams::fgkSampleMin, min = AliEMCALGeoParams::fgkSampleMax; // min and max sample values
474       int nsamples = 0;
475
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++) {
480           sample = sig[i];
481
482           // check if it's a min or max value
483           if (sample < min) min = sample;
484           if (sample > max) max = sample;
485
486         } // loop over samples in bunch
487       } // loop over bunches
488
489       if (nsamples > 0) { // this check is needed for when we have zero-supp. on, but not sparse readout
490
491       gain = -1; // init to not valid value
492       //If we're here then we're done with this tower
493       if ( in->IsLowGain() ) {
494         gain = 0;
495       }
496       else if ( in->IsHighGain() ) {
497         gain = 1;
498       }
499       else if ( in->IsLEDMonData() ) {
500         gain = in->GetRow(); // gain coded in (in RCU/Altro mapping) as Row info for LED refs..
501       }
502
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
505       //Debug
506       if (arrayPos < 0 || arrayPos >= fModules) {
507         printf("AliCaloCalibSignal::ProcessEvent = Oh no: arrayPos = %i.\n", arrayPos); 
508         return kFALSE;
509       }
510
511       if ( in->IsHighGain() || in->IsLowGain() ) { // regular tower
512         // get tower number for AmpVal array
513         iTowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow()); 
514
515         if (gain == 0) {
516           // fill amplitude into the array         
517           iAmpValLowGain[iTowerNum] = max - min;
518           nLowChan++;
519         } 
520         else if (gain==1) {//fill the high gain ones
521           // fill amplitude into the array
522           iAmpValHighGain[iTowerNum] = max - min;
523           nHighChan++;
524         }//end if gain
525       } // regular tower
526       else if ( in->IsLEDMonData() ) { // LED ref.; 
527         // strip # is coded is 'column' in the channel maps 
528         iRefNum = GetRefNum(arrayPos, in->GetColumn(), gain); 
529         iLEDAmpVal[iRefNum] = max - min;
530       } // end of LED ref
531
532       } // nsamples>0 check, some data found for this channel; not only trailer/header      
533     } // end while over channel 
534    
535   }//end while over DDL's, of input stream
536   
537   in->Reset(); // just in case the next customer forgets to check if the stream was reset..
538
539   // now check if it was a led event, only use high gain (that should be sufficient)
540
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; 
545   }
546
547   if (fReqFractionAboveAmp) {
548     bool ok = false;
549     if (nHighChan > 0) { 
550       ok = CheckFractionAboveAmp(iAmpValHighGain, checkResultArray); 
551     }
552     if (!ok) return false; // skip event
553   }
554
555   fNAcceptedEvents++; // one more event accepted
556
557   if (fStartTime == 0) { // if start-timestamp wasn't set,we'll pick it up from the first event we encounter
558     fStartTime = Timestamp;
559   }
560
561   fHour = (Timestamp - fStartTime)/(double)fgkNumSecInHr;
562   if (fLatestHour < fHour) {
563     fLatestHour = fHour; 
564   }
565   
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++){
571         
572         iTowerNum = GetTowerNum(i, j, k); 
573
574         if(iAmpValHighGain[iTowerNum]) {
575           fAmp = iAmpValHighGain[iTowerNum];
576           fChannelNum = GetChannelNum(i,j,k,1);
577           fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValHighGain[iTowerNum]);
578           fNHighGain[iTowerNum]++;
579         }
580         if(iAmpValLowGain[iTowerNum]) {
581           fAmp = iAmpValLowGain[iTowerNum];
582           fChannelNum = GetChannelNum(i,j,k,0);
583           fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValLowGain[iTowerNum]);
584           fNLowGain[iTowerNum]++;
585         }
586       } // rows
587       } // column passed check
588     } // columns
589
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;
594     // which leads to
595     // iColFirst = (iSM%2==0) ? iStrip*2 : (AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iStrip)*2;
596
597     // also LED refs
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 (iLEDAmpVal[fRefNum]) {
604           fAmp = iLEDAmpVal[fRefNum];
605           fTreeLEDAmpVsTime->Fill();//fRefNum,fHour,fAmp);
606           fNRef[fRefNum]++;
607         }
608       } // gain
609       } // at least one column in strip passed check 
610     }
611
612   } // modules
613   
614   return kTRUE;
615 }
616
617 //_____________________________________________________________________
618 Bool_t AliCaloCalibSignal::Save(TString fileName)
619 {
620   //Saves all the TTrees to the designated file
621   
622   TFile destFile(fileName, "recreate");
623   
624   if (destFile.IsZombie()) {
625     return kFALSE;
626   }
627   
628   destFile.cd();
629
630   // save the trees
631   fTreeAmpVsTime->Write();
632   fTreeLEDAmpVsTime->Write();
633   if (fUseAverage) { 
634     Analyze(); // get the latest and greatest averages
635     fTreeAvgAmpVsTime->Write();
636     fTreeLEDAvgAmpVsTime->Write();
637   }
638
639   destFile.Close();
640   
641   return kTRUE;
642 }
643
644 //_____________________________________________________________________
645 Bool_t AliCaloCalibSignal::Analyze()
646 {
647   // Fill the tree holding the average values
648   if (!fUseAverage) { return kFALSE; }
649
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();
654   }
655
656   //0: setup variables for the TProfile plots that we'll use to do the averages
657   int numProfBins = 0;
658   double timeMin = 0;
659   double timeMax = 0;
660   if (fSecInAverage > 0) {
661     numProfBins = (int)( (fLatestHour*fgkNumSecInHr)/fSecInAverage + 1 ); // round-off
662   }
663   numProfBins += 2; // add extra buffer : first and last
664   double binSize = 1.0*fSecInAverage / fgkNumSecInHr;
665   timeMin = - binSize;
666   timeMax = timeMin + numProfBins*binSize;
667
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));
671
672   char name[200]; // for profile id and title
673   int iTowerNum = 0;
674
675   for (int i = 0; i<fModules; i++) {
676     for (int ic=0; ic<fColumns; ic++){
677       for (int ir=0; ir<fRows; ir++) {
678
679         iTowerNum = GetTowerNum(i, ic, ir);
680         // high gain
681         if (fNHighGain[iTowerNum] > 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");
685         }
686
687         // same for low gain
688         if (fNLowGain[iTowerNum] > 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");
692         }
693
694       } // rows
695     } // columns
696   } // modules
697
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);
703
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);
709     }
710   }
711
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);
717
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
729         } // loop over bins
730       } // some entries for this profile
731     } // profile exists  
732   } // loop over all possible channels
733
734
735   // and finally, go through same exercise for LED also.. 
736
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));
740
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");
748         } 
749       }// gain
750     } 
751   } // modules
752
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);
758
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);
764     }
765   }
766
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);
772
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
784         } // loop over bins
785       } // some entries for this profile
786     } // profile exists  
787   } // loop over all possible channels
788
789   // OK, we're done..
790
791   return kTRUE;
792 }
793
794 //_____________________________________________________________________
795 Bool_t AliCaloCalibSignal::DecodeChannelNum(const int chanId, 
796                                             int *imod, int *icol, int *irow, int *igain) const  
797 { // return the module, column, row, and gain for a given channel number
798   *igain = chanId/(fModules*fColumns*fRows);
799   *imod = (chanId/(fColumns*fRows)) % fModules;
800   *icol = (chanId/fRows) % fColumns;
801   *irow = chanId % fRows;
802   return kTRUE;
803
804
805 //_____________________________________________________________________
806 Bool_t AliCaloCalibSignal::DecodeRefNum(const int refId, 
807                                         int *imod, int *istripMod, int *igain) const 
808 { // return the module, stripModule, and gain for a given reference number
809   *igain = refId/(fModules*fLEDRefs);
810   *imod = (refId/(fLEDRefs)) % fModules;
811   *istripMod = refId % fLEDRefs;
812   return kTRUE;
813