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