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