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