AliRoot Core  ee782a0 (ee782a0)
AliCaloCalibSignal.cxx
Go to the documentation of this file.
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 
16 // --- Standard library ---
17 #include <string>
18 #include <sstream>
19 #include <fstream>
20 
21 // --- ROOT system ---
22 #include "TProfile.h"
23 #include "TFile.h"
24 
25 // --- AliRoot header files ---
26 #include "AliRawReader.h"
27 #include "AliCaloRawStreamV3.h"
28 #include "AliDAQ.h"
29 
30 // --- Calo header files ---
31 #include "AliCaloConstants.h"
32 #include "AliCaloBunchInfo.h"
33 #include "AliCaloFitResults.h"
34 #include "AliCaloRawAnalyzer.h"
36 
37 #include "AliCaloCalibSignal.h"
38 
40 ClassImp(AliCaloCalibSignal) ;
42 
43 using namespace std;
44 
45 // variables for TTree filling; not sure if they should be static or not
46 static int fChannelNum = 0;
47 static int fRefNum = 0;
48 static double fAmp = 0;
49 static double fAvgAmp = 0;
50 static double fRMS = 0;
51 
56 //_____________________________________________________________________
58  TObject(),
59  fDetType(kNone),
60  fColumns(0),
61  fRows(0),
62  fLEDRefs(0),
63  fModules(0),
64  fCaloString(),
65  fMapping(NULL),
66  fFittingAlgorithm(0),
67  fRawAnalyzer(0),
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  fDownscale(10),
80  fNEvents(0),
81  fNAcceptedEvents(0),
82  fTreeAmpVsTime(NULL),
83  fTreeAvgAmpVsTime(NULL),
84  fTreeLEDAmpVsTime(NULL),
85  fTreeLEDAvgAmpVsTime(NULL)
86 {
87  // First we set the detector-type related constants.
88  if (detectorType == kPhos)
89  {
94  fCaloString = "PHOS";
95  }
96  else
97  {
98  // We'll just trust the enum to keep everything in line, so that if detectorType
99  // isn't kPhos then it is kEmCal. Note, however, that this is not necessarily the
100  // case, if someone intentionally gives another number
104  fModules = 20; // AliEMCALGeoParams::fgkEMCALModules is set to a higher number due to simulation constraints...
105  fCaloString = "EMCAL";
106  }
107 
108  fDetType = detectorType;
110  ResetInfo(); // trees and counters
111 }
112 
115 //_____________________________________________________________________
117 {
118  DeleteTrees();
119 }
120 
123 //_____________________________________________________________________
125 {
126  if (fTreeAmpVsTime) delete fTreeAmpVsTime;
130 
131  // and reset pointers
132  fTreeAmpVsTime = NULL;
133  fTreeAvgAmpVsTime = NULL;
134  fTreeLEDAmpVsTime = NULL;
135  fTreeLEDAvgAmpVsTime = NULL;
136 
137  return;
138 }
139 
140 // copy ctor
141 //_____________________________________________________________________
142 //AliCaloCalibSignal::AliCaloCalibSignal(const AliCaloCalibSignal &sig) :
143 // TObject(sig),
144 // fDetType(sig.GetDetectorType()),
145 // fColumns(sig.GetColumns()),
146 // fRows(sig.GetRows()),
147 // fLEDRefs(sig.GetLEDRefs()),
148 // fModules(sig.GetModules()),
149 // fCaloString(sig.GetCaloString()),
150 // fMapping(), //! note that we are not copying the map info
151 // fRunNumber(sig.GetRunNumber()),
152 // fStartTime(sig.GetStartTime()),
153 // fAmpCut(sig.GetAmpCut()),
154 // fReqFractionAboveAmpCutVal(sig.GetReqFractionAboveAmpCutVal()),
155 // fReqFractionAboveAmp(sig.GetReqFractionAboveAmp()),
156 // fAmpCutLEDRef(sig.GetAmpCutLEDRef()),
157 // fReqLEDRefAboveAmpCutVal(sig.GetReqLEDRefAboveAmpCutVal()),
158 // fHour(sig.GetHour()),
159 // fLatestHour(sig.GetLatestHour()),
160 // fUseAverage(sig.GetUseAverage()),
161 // fSecInAverage(sig.GetSecInAverage()),
162 // fDownscale(sig.GetDownscale()),
163 // fNEvents(sig.GetNEvents()),
164 // fNAcceptedEvents(sig.GetNAcceptedEvents()),
165 // fTreeAmpVsTime(),
166 // fTreeAvgAmpVsTime(),
167 // fTreeLEDAmpVsTime(),
168 // fTreeLEDAvgAmpVsTime()
169 //{
170 // // also the TTree contents
171 // AddInfo(&sig);
172 // for (Int_t i = 0; i<fgkMaxTowers; i++) {
173 // fNHighGain[i] = sig.fNHighGain[i];
174 // fNLowGain[i] = sig.fNLowGain[i];
175 // }
176 // for (Int_t i = 0; i<(2*fgkMaxRefs); i++) {
177 // fNRef[i] = sig.fNRef[i];
178 // }
179 //
180 //
181 //}
182 //
183 // assignment operator; use copy ctor to make life easy..
184 //_____________________________________________________________________
185 //AliCaloCalibSignal& AliCaloCalibSignal::operator = (const AliCaloCalibSignal &source)
186 //{
187 // // assignment operator; use copy ctor
188 // if (&source == this) return *this;
189 //
190 // new (this) AliCaloCalibSignal(source);
191 // return *this;
192 //}
193 
196 //_____________________________________________________________________
198 {
199  // first, regular version
200  fTreeAmpVsTime = new TTree("fTreeAmpVsTime","Amplitude vs. Time Tree Variables");
201 
202  fTreeAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I");
203  fTreeAmpVsTime->Branch("fHour", &fHour, "fHour/D");
204  fTreeAmpVsTime->Branch("fAmp", &fAmp, "fAmp/D");
205 
206  // then, average version
207  fTreeAvgAmpVsTime = new TTree("fTreeAvgAmpVsTime","Average Amplitude vs. Time Tree Variables");
208 
209  fTreeAvgAmpVsTime->Branch("fChannelNum", &fChannelNum, "fChannelNum/I");
210  fTreeAvgAmpVsTime->Branch("fHour", &fHour, "fHour/D");
211  fTreeAvgAmpVsTime->Branch("fAvgAmp", &fAvgAmp, "fAvgAmp/D");
212  fTreeAvgAmpVsTime->Branch("fRMS", &fRMS, "fRMS/D");
213 
214  // then same for LED..
215  fTreeLEDAmpVsTime = new TTree("fTreeLEDAmpVsTime","LED Amplitude vs. Time Tree Variables");
216  fTreeLEDAmpVsTime->Branch("fRefNum", &fRefNum, "fRefNum/I");
217  fTreeLEDAmpVsTime->Branch("fHour", &fHour, "fHour/D");
218  fTreeLEDAmpVsTime->Branch("fAmp", &fAmp, "fAmp/D");
219 
220  fTreeLEDAvgAmpVsTime = new TTree("fTreeLEDAvgAmpVsTime","Average LED Amplitude vs. Time Tree Variables");
221  fTreeLEDAvgAmpVsTime->Branch("fRefNum", &fRefNum, "fRefNum/I");
222  fTreeLEDAvgAmpVsTime->Branch("fHour", &fHour, "fHour/D");
223  fTreeLEDAvgAmpVsTime->Branch("fAvgAmp", &fAvgAmp, "fAvgAmp/D");
224  fTreeLEDAvgAmpVsTime->Branch("fRMS", &fRMS, "fRMS/D");
225 
226  return;
227 }
228 
231 //_____________________________________________________________________
233 {
234  Zero(); // set all counters to 0
235  DeleteTrees(); // delete previous stuff
236  CreateTrees(); // and create some new ones
237 
238  return;
239 }
240 
243 //_____________________________________________________________________
245 {
246  fHour = 0;
247  fLatestHour = 0;
248  fNEvents = 0;
249  fNAcceptedEvents = 0;
250 
251  // Set the number of points for each tower: Amp vs. Time
252  memset(fNHighGain, 0, sizeof(fNHighGain));
253  memset(fNLowGain, 0, sizeof(fNLowGain));
254  // and LED reference
255  memset(fNRef, 0, sizeof(fNRef));
256 
257  return;
258 }
259 
263 //_____________________________________________________________________
264 Bool_t AliCaloCalibSignal::CheckFractionAboveAmp(const int *iAmpVal,
265  int resultArray[]) const
266 {
267  Bool_t returnCode = false;
268 
269  int iTowerNum = 0;
270  double fraction = 0;
271  for (int i = 0; i<fModules; i++)
272  {
273  for (int j = 0; j<fColumns; j++)
274  {
275  Int_t nRows = fRows;
276  if (i==10 || i==11 || i==18 || i==19)
277  nRows=fRows/3;
278  int nAbove = 0;
279  for (int k = 0; k<fRows; k++)
280  {
281  iTowerNum = GetTowerNum(i,j,k);
282  if (iAmpVal[iTowerNum] > fAmpCut)
283  {
284  nAbove++;
285  }
286  }
287 
288  resultArray[i*fColumns +j] = 0; // init. to denied
289 
290  if (nAbove > 0)
291  {
292  fraction = (1.0*nAbove) / nRows;
293  /*
294  printf("DS mod %d col %d nAbove %d fraction %3.2f\n",
295  i, j, nAbove, fraction);
296  */
297  if (fraction > fReqFractionAboveAmpCutVal)
298  {
299  resultArray[i*fColumns + j] = nAbove;
300  returnCode = true;
301  }
302  }
303  }
304  } // modules loop
305 
306  return returnCode;
307 }
308 
312 //_____________________________________________________________________
313 Bool_t AliCaloCalibSignal::CheckLEDRefAboveAmp(const int *iAmpVal,
314  int resultArray[]) const
315 {
316  Bool_t returnCode = false;
317 
318  int iRefNum = 0;
319  int gain = 1; // look at high gain; this should be rather saturated usually..
320  for (int i = 0; i<fModules; i++)
321  {
322  for (int j = 0; j<fLEDRefs; j++)
323  {
324  iRefNum = GetRefNum(i, j, gain);
325  if (iAmpVal[iRefNum] > fAmpCutLEDRef)
326  {
327  resultArray[i*fLEDRefs +j] = 1; // enough signal
328  returnCode = true;
329  }
330  else
331  {
332  resultArray[i*fLEDRefs +j] = 0; // not enough signal
333  }
334 
335  /*
336  printf("DS mod %d LEDRef %d ampVal %d\n",
337  i, j, iAmpVal[iRefNum]);
338  */
339  } // LEDRefs
340  } // modules loop
341 
342  return returnCode;
343 }
344 
349 //_____________________________________________________________________
350 void AliCaloCalibSignal::SetParametersFromFile(const char *parameterFile)
351 {
352  static const string delimitor("::");
353 
354  // open, check input file
355  ifstream in( parameterFile );
356  if( !in )
357  {
358  printf("in AliCaloCalibSignal::SetParametersFromFile - Using default/run_time parameters.\n");
359  return;
360  }
361 
362  // Note: this method is a bit more complicated than it really has to be
363  // - allowing for multiple entries per line, arbitrary order of the
364  // different variables etc. But I wanted to try and do this in as
365  // correct a C++ way as I could (as an exercise).
366 
367  // read in
368  char readline[1024];
369  while ((in.rdstate() & ios::failbit) == 0 )
370  {
371  // Read into the raw char array and then construct a string
372  // to do the searching
373  in.getline(readline, 1024);
374  istringstream s(readline);
375 
376  while ( ( s.rdstate() & ios::failbit ) == 0 )
377  {
378  string keyValue;
379  s >> keyValue;
380 
381  // check stream status
382  if( ( s.rdstate() & ios::failbit ) == ios::failbit ) break;
383 
384  // skip rest of line if comments found
385  if( keyValue.substr( 0, 2 ) == "//" ) break;
386 
387  // look for "::" in keyValue pair
388  size_t position = keyValue.find( delimitor );
389  if( position == string::npos )
390  {
391  printf("wrong format for key::value pair: %s\n", keyValue.c_str());
392  }
393 
394  // split keyValue pair
395  string key( keyValue.substr( 0, position ) );
396  string value( keyValue.substr( position+delimitor.size(),
397  keyValue.size()-delimitor.size() ) );
398 
399  // check value does not contain a new delimitor
400  if( value.find( delimitor ) != string::npos )
401  {
402  printf("wrong format for key::value pair: %s\n", keyValue.c_str());
403  }
404 
405  // debug: check key value pair
406  // printf("AliCaloCalibSignal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
407 
408  // if the key matches with something we expect, we assign the new value
409  if ( (key == "fAmpCut") || (key == "fReqFractionAboveAmpCutVal") ||
410  (key == "fAmpCutLEDRef") || (key == "fSecInAverage") ||
411  (key == "fFittingAlgorithm") || (key == "fDownscale") )
412  {
413  istringstream iss(value);
414  printf("AliCaloCalibSignal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
415 
416  if (key == "fAmpCut") {
417  iss >> fAmpCut;
418  }
419  else if (key == "fReqFractionAboveAmpCutVal") {
421  }
422  else if (key == "fAmpCutLEDRef") {
423  iss >> fAmpCutLEDRef;
424  }
425  else if (key == "fSecInAverage") {
426  iss >> fSecInAverage;
427  }
428  else if (key == "fFittingAlgorithm") {
429  iss >> fFittingAlgorithm;
430  SetFittingAlgorithm( fFittingAlgorithm );
431  }
432  else if (key == "fDownscale") {
433  iss >> fDownscale;
434  }
435  } // some match found/expected
436 
437  }
438  }
439 
440  in.close();
441  return;
442 }
443 
446 //_____________________________________________________________________
447 void AliCaloCalibSignal::WriteParametersToFile(const char *parameterFile)
448 {
449  static const string delimitor("::");
450 
451  ofstream out( parameterFile );
452  out << "// " << parameterFile << endl;
453  out << "fAmpCut" << "::" << fAmpCut << endl;
454  out << "fReqFractionAboveAmpCutVal" << "::" << fReqFractionAboveAmpCutVal << endl;
455  out << "fAmpCutLEDRef" << "::" << fAmpCutLEDRef << endl;
456  out << "fSecInAverage" << "::" << fSecInAverage << endl;
457  out << "fFittingAlgorithm" << "::" << fFittingAlgorithm << endl;
458  out << "fDownscale" << "::" << fDownscale << endl;
459 
460  out.close();
461  return;
462 }
463 
468 //_____________________________________________________________________
470 {
471  fFittingAlgorithm = fitAlgo;
472 
473  delete fRawAnalyzer; // delete doesn't do anything if the pointer is 0x0
474 
476  fRawAnalyzer->SetIsZeroSuppressed(true); // TMP - should use stream->IsZeroSuppressed(), or altro cfg registers later
477 }
478 
483 //_____________________________________________________________________
485 {
486  // note/FIXME: we are not yet adding correctly the info for fN{HighGain,LowGain,Ref}
487  // here - but consider this a feature for now (20080905):
488  // we'll do Analyze() unless entries were found for a tower in this original object.
489 
490  // add info from sig's TTrees to ours..
491  TTree *sigAmp = sig->GetTreeAmpVsTime();
492  TTree *sigAvgAmp = sig->GetTreeAvgAmpVsTime();
493 
494  // we could try some merging via TList or what also as a more elegant approach
495  // but I wanted with the stupid/simple and hopefully safe approach of looping
496  // over what we want to add..
497 
498  // associate variables for sigAmp and sigAvgAmp:
499  sigAmp->SetBranchAddress("fChannelNum",&fChannelNum);
500  sigAmp->SetBranchAddress("fHour",&fHour);
501  sigAmp->SetBranchAddress("fAmp",&fAmp);
502 
503  // loop over the trees.. note that since we use the same variables we should not need
504  // to do any assignments between the getting and filling
505  for (int i=0; i<sigAmp->GetEntries(); i++)
506  {
507  sigAmp->GetEntry(i);
508  fTreeAmpVsTime->Fill();
509  }
510 
511  sigAvgAmp->SetBranchAddress("fChannelNum",&fChannelNum);
512  sigAvgAmp->SetBranchAddress("fHour",&fHour);
513  sigAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
514  sigAvgAmp->SetBranchAddress("fRMS",&fRMS);
515 
516  for (int i=0; i<sigAvgAmp->GetEntries(); i++)
517  {
518  sigAvgAmp->GetEntry(i);
519  fTreeAvgAmpVsTime->Fill();
520  }
521 
522  // also LED..
523  TTree *sigLEDAmp = sig->GetTreeLEDAmpVsTime();
524  TTree *sigLEDAvgAmp = sig->GetTreeLEDAvgAmpVsTime();
525 
526  // associate variables for sigAmp and sigAvgAmp:
527  sigLEDAmp->SetBranchAddress("fRefNum",&fRefNum);
528  sigLEDAmp->SetBranchAddress("fHour",&fHour);
529  sigLEDAmp->SetBranchAddress("fAmp",&fAmp);
530 
531  // loop over the trees.. note that since we use the same variables we should not need
532  // to do any assignments between the getting and filling
533  for (int i=0; i<sigLEDAmp->GetEntries(); i++)
534  {
535  sigLEDAmp->GetEntry(i);
536  fTreeLEDAmpVsTime->Fill();
537  }
538 
539  sigLEDAvgAmp->SetBranchAddress("fRefNum",&fRefNum);
540  sigLEDAvgAmp->SetBranchAddress("fHour",&fHour);
541  sigLEDAvgAmp->SetBranchAddress("fAvgAmp",&fAvgAmp);
542  sigLEDAvgAmp->SetBranchAddress("fRMS",&fRMS);
543 
544  for (int i=0; i<sigLEDAvgAmp->GetEntries(); i++)
545  {
546  sigLEDAvgAmp->GetEntry(i);
547  fTreeLEDAvgAmpVsTime->Fill();
548  }
549 
550  // We should also copy other pieces of info: counters and parameters
551  // (not number of columns and rows etc which should be the same)
552  // note that I just assign them here rather than Add them, but we
553  // normally just Add (e.g. in Preprocessor) one object so this should be fine.
554  fRunNumber = sig->GetRunNumber();
555  fStartTime = sig->GetStartTime();
556  fAmpCut = sig->GetAmpCut();
561  fHour = sig->GetHour();
562  fLatestHour = sig->GetLatestHour();
563  fUseAverage = sig->GetUseAverage();
565  fDownscale = sig->GetDownscale();
566  fNEvents = sig->GetNEvents();
568 
569  return kTRUE;//We hopefully succesfully added info from the supplied object
570 }
571 
576 //_____________________________________________________________________
577 Bool_t AliCaloCalibSignal::ProcessEvent(AliRawReader *rawReader)
578 {
579  // if fMapping is NULL the rawstream will crate its own mapping
580  AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
581  if (fDetType == kEmCal)
582  rawReader->Select("EMCAL",0,AliDAQ::GetFirstSTUDDL()-1) ; //select EMCAL DDL range
583 
584  return ProcessEvent( &rawStream, rawReader->GetTimestamp() );
585 }
586 
592 //_____________________________________________________________________
593 Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStreamV3 *in, UInt_t timestamp)
594 {
595  if (!in) return kFALSE; //Return right away if there's a null pointer
596 
597  fNEvents++; // one more event
598 
599  if ( (fNEvents%fDownscale)!=0 ) return kFALSE; // mechanism to skip some of the input events, if we want
600 
601  // use maximum numbers to set array sizes
602  int iAmpValHighGain[fgkMaxTowers];
603  int iAmpValLowGain[fgkMaxTowers];
604  memset(iAmpValHighGain, 0, sizeof(iAmpValHighGain));
605  memset(iAmpValLowGain, 0, sizeof(iAmpValLowGain));
606 
607  // also for LED reference
608  int iLEDAmpVal[fgkMaxRefs * 2]; // factor 2 is for the two gain values
609  memset(iLEDAmpVal, 0, sizeof(iLEDAmpVal));
610 
611  int gain = 0; // high or low gain
612 
613  // Number of Low and High gain, and LED Ref, channels for this event:
614  int nLowChan = 0;
615  int nHighChan = 0;
616  int nLEDRefChan = 0;
617 
618  int iTowerNum = 0; // array index for regular towers
619  int iRefNum = 0; // array index for LED references
620 
621  // loop first to get the fraction of channels with amplitudes above cut
622 
623  while (in->NextDDL())
624  {
625  while (in->NextChannel())
626  {
627  vector<AliCaloBunchInfo> bunchlist;
628  while (in->NextBunch())
629  {
630  bunchlist.push_back( AliCaloBunchInfo(in->GetStartTimeBin(), in->GetBunchLength(), in->GetSignals() ) );
631  }
632 
633  if (bunchlist.size() == 0) continue;
634 
635  gain = -1; // init to not valid value
636  //If we're here then we're done with this tower
637  if ( in->IsLowGain() )
638  {
639  gain = 0;
640  }
641  else if ( in->IsHighGain() )
642  {
643  gain = 1;
644  }
645  else if ( in->IsLEDMonData() )
646  {
647  gain = in->GetRow(); // gain coded in (in RCU/Altro mapping) as Row info for LED refs..
648  }
649  else { continue ; } // don't try to fit TRU..
650 
651  // it should be enough to check the SuperModule info for each DDL really, but let's keep it here for now
652  int arrayPos = in->GetModule(); //The modules are numbered starting from 0
653  //Debug
654  if (arrayPos < 0 || arrayPos >= fModules)
655  {
656  printf("AliCaloCalibSignal::ProcessEvent = Oh no: arrayPos = %i.\n", arrayPos);
657  return kFALSE;
658  }
659 
660  AliCaloFitResults res = fRawAnalyzer->Evaluate( bunchlist, in->GetAltroCFG1(), in->GetAltroCFG2());
661  if ( in->IsHighGain() || in->IsLowGain() )
662  { // regular tower
663  // get tower number for AmpVal array
664  iTowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow());
665 
666  if (gain == 0)
667  {
668  // fill amplitude into the array
669  iAmpValLowGain[iTowerNum] = (int) res.GetAmp();
670  nLowChan++;
671  }
672  else if (gain==1)
673  {//fill the high gain ones
674  // fill amplitude into the array
675  iAmpValHighGain[iTowerNum] = (int) res.GetAmp();
676  nHighChan++;
677  }//end if gain
678  } // regular tower
679  else if ( in->IsLEDMonData() )
680  { // LED ref.;
681  // strip # is coded is 'column' in the channel maps
682  iRefNum = GetRefNum(arrayPos, in->GetColumn(), gain);
683  iLEDAmpVal[iRefNum] = (int) res.GetAmp();
684  nLEDRefChan++;
685  } // end of LED ref
686 
687  } // end while over channel
688 
689  }//end while over DDL's, of input stream
690 
691  in->Reset(); // just in case the next customer forgets to check if the stream was reset..
692 
693  // now check if it was an LED event, using the LED Reference info per strip
694 
695  // by default all columns are accepted (init check to > 0)
698  {
699  checkResultArray[ia] = 1;
700  }
701 
703  {
704  bool ok = false;
705  if (nHighChan > 0)
706  {
707  ok = CheckFractionAboveAmp(iAmpValHighGain, checkResultArray);
708  }
709 
710  if (!ok) return false; // skip event
711  }
712 
713  // by default all columns are accepted (init check to > 0)
716  {
717  checkResultArrayLEDRef[ia] = 1;
718  }
719 
721  {
722  bool ok = false;
723  if (nLEDRefChan > 0)
724  {
725  ok = CheckLEDRefAboveAmp(iLEDAmpVal, checkResultArrayLEDRef);
726  }
727 
728  if (!ok) return false; // skip event
729  }
730 
731  fNAcceptedEvents++; // one more event accepted
732 
733  if (fStartTime == 0)
734  { // if start-timestamp wasn't set,we'll pick it up from the first event we encounter
735  fStartTime = timestamp;
736  }
737 
738  fHour = (timestamp - fStartTime)/(double)fgkNumSecInHr;
739  if (fLatestHour < fHour)
740  {
741  fLatestHour = fHour;
742  }
743 
744  // it is a led event, now fill TTree
745  // We also do the activity check for LEDRefs/Strips, but need to translate between column
746  // and strip indices for that; based on these relations:
747  // iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
748  // iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2;
749  // which leads to
750  // iColFirst = (iSM%2==0) ? iStrip*2 : (AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iStrip)*2;
751 
752  for(int i=0; i<fModules; i++)
753  {
754  for(int j=0; j<fColumns; j++)
755  {
756  int iStrip = (i%2==0) ? j/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - j/2;
757  if (checkResultArray[i*fColumns + j]>0 && checkResultArrayLEDRef[i*fLEDRefs + iStrip]>0)
758  {
759  // column passed check
760  for(int k=0; k<fRows; k++)
761  {
762  iTowerNum = GetTowerNum(i, j, k);
763 
764  if(iAmpValHighGain[iTowerNum])
765  {
766  fAmp = iAmpValHighGain[iTowerNum];
767  fChannelNum = GetChannelNum(i,j,k,1);
768  fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValHighGain[iTowerNum]);
769  fNHighGain[iTowerNum]++;
770  }
771 
772  if(iAmpValLowGain[iTowerNum])
773  {
774  fAmp = iAmpValLowGain[iTowerNum];
775  fChannelNum = GetChannelNum(i,j,k,0);
776  fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValLowGain[iTowerNum]);
777  fNLowGain[iTowerNum]++;
778  }
779  } // rows
780  } // column passed check, and LED Ref for strip passed check (if any)
781  } // columns
782 
783  // also LED refs
784  for(int j=0; j<fLEDRefs; j++)
785  {
786  int iColFirst = (i%2==0) ? j*2 : (AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - j)*2; //CHECKME!!!
787  if ( ((checkResultArray[i*fColumns + iColFirst]>0) || (checkResultArray[i*fColumns + iColFirst + 1]>0)) && // at least one column in strip passed check
788  (checkResultArrayLEDRef[i*fLEDRefs + j]>0) )
789  { // and LED Ref passed checks
790  for (gain=0; gain<2; gain++)
791  {
792  fRefNum = GetRefNum(i, j, gain);
793  if (iLEDAmpVal[fRefNum])
794  {
795  fAmp = iLEDAmpVal[fRefNum];
796  fTreeLEDAmpVsTime->Fill();//fRefNum,fHour,fAmp);
797  fNRef[fRefNum]++;
798  }
799  } // gain
800  } // at least one column in strip passed check, and LED Ref passed check (if any)
801  }
802 
803  } // modules
804 
805  return kTRUE;
806 }
807 
812 //_____________________________________________________________________
814 {
815  TFile destFile(fileName, "recreate");
816 
817  if (destFile.IsZombie())
818  return kFALSE;
819 
820  destFile.cd();
821 
822  // save the trees
823  fTreeAmpVsTime->Write();
824  fTreeLEDAmpVsTime->Write();
825 
826  if (fUseAverage)
827  {
828  Analyze(); // get the latest and greatest averages
829  fTreeAvgAmpVsTime->Write();
830  fTreeLEDAvgAmpVsTime->Write();
831  }
832 
833  destFile.Close();
834 
835  return kTRUE;
836 }
837 
841 //_____________________________________________________________________
843 {
844  if (!fUseAverage) { return kFALSE; }
845 
846  // Reset the average TTree if Analyze has already been called earlier,
847  // meaning that the TTree could have been partially filled
848  if (fTreeAvgAmpVsTime->GetEntries() > 0)
849  fTreeAvgAmpVsTime->Reset();
850 
851  //0: setup variables for the TProfile plots that we'll use to do the averages
852  int numProfBins = 0;
853  double timeMin = 0;
854  double timeMax = 0;
855  if (fSecInAverage > 0)
856  numProfBins = (int)( (fLatestHour*fgkNumSecInHr)/fSecInAverage + 1 ); // round-off
857 
858  numProfBins += 2; // add extra buffer : first and last
859  double binSize = 1.0*fSecInAverage / fgkNumSecInHr;
860  timeMin = - binSize;
861  timeMax = timeMin + numProfBins*binSize;
862 
863  //1: set up TProfiles for the towers that had data
864  TProfile * profile[fgkMaxTowers*2]; // *2 is since we include both high and low gains
865  memset(profile, 0, sizeof(profile));
866  const Int_t buffersize = 200;
867  char name[buffersize]; // for profile id and title
868  int iTowerNum = 0;
869 
870  for (int i = 0; i<fModules; i++)
871  {
872  for (int ic=0; ic<fColumns; ic++)
873  {
874  for (int ir=0; ir<fRows; ir++)
875  {
876  iTowerNum = GetTowerNum(i, ic, ir);
877  // high gain
878  if (fNHighGain[iTowerNum] > 0)
879  {
880  fChannelNum = GetChannelNum(i, ic, ir, 1);
881  snprintf(name,buffersize,"profileChan%d", fChannelNum);
882  profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
883  }
884 
885  // same for low gain
886  if (fNLowGain[iTowerNum] > 0)
887  {
888  fChannelNum = GetChannelNum(i, ic, ir, 0);
889  snprintf(name,buffersize,"profileChan%d", fChannelNum);
890  profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
891  }
892 
893  } // rows
894  } // columns
895  } // modules
896 
897  //2: fill profiles by looping over tree
898  // Set addresses for tree-readback also
899  fTreeAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
900  fTreeAmpVsTime->SetBranchAddress("fHour", &fHour);
901  fTreeAmpVsTime->SetBranchAddress("fAmp", &fAmp);
902 
903  for (int ient=0; ient<fTreeAmpVsTime->GetEntries(); ient++)
904  {
905  fTreeAmpVsTime->GetEntry(ient);
906  if (profile[fChannelNum])
907  {
908  // profile should always have been created above, for active channels
909  profile[fChannelNum]->Fill(fHour, fAmp);
910  }
911  }
912 
913  // re-associating the branch addresses here seems to be needed for OK 'average' storage
914  fTreeAvgAmpVsTime->SetBranchAddress("fChannelNum", &fChannelNum);
915  fTreeAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
916  fTreeAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
917  fTreeAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);
918 
919  //3: fill avg tree by looping over the profiles
921  {
922  if (profile[fChannelNum]) { // profile was created
923  if (profile[fChannelNum]->GetEntries() > 0)
924  {
925  // profile had some entries
926  for(int it=0; it<numProfBins; it++)
927  {
928  if (profile[fChannelNum]->GetBinEntries(it+1) > 0)
929  {
930  fAvgAmp = profile[fChannelNum]->GetBinContent(it+1);
931  fHour = profile[fChannelNum]->GetBinCenter(it+1);
932  fRMS = profile[fChannelNum]->GetBinError(it+1);
933  fTreeAvgAmpVsTime->Fill();
934  } // some entries for this bin
935  } // loop over bins
936  } // some entries for this profile
937  } // profile exists
938  } // loop over all possible channels
939 
940 
941  // and finally, go through same exercise for LED also..
942 
943  //1: set up TProfiles for the towers that had data
944  TProfile * profileLED[fgkMaxRefs*2]; // *2 is since we include both high and low gains
945  memset(profileLED, 0, sizeof(profileLED));
946 
947  for (int i = 0; i<fModules; i++)
948  {
949  for(int j=0; j<fLEDRefs; j++)
950  {
951  for (int gain=0; gain<2; gain++)
952  {
953  fRefNum = GetRefNum(i, j, gain);
954  if (fNRef[fRefNum] > 0)
955  {
956  snprintf(name, buffersize, "profileLEDRef%d", fRefNum);
957  profileLED[fRefNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
958  }
959  }// gain
960  }
961  } // modules
962 
963  //2: fill profiles by looping over tree
964  // Set addresses for tree-readback also
965  fTreeLEDAmpVsTime->SetBranchAddress("fRefNum", &fRefNum);
966  fTreeLEDAmpVsTime->SetBranchAddress("fHour", &fHour);
967  fTreeLEDAmpVsTime->SetBranchAddress("fAmp", &fAmp);
968 
969  for (int ient=0; ient<fTreeLEDAmpVsTime->GetEntries(); ient++)
970  {
971  fTreeLEDAmpVsTime->GetEntry(ient);
972  if (profileLED[fRefNum])
973  {
974  // profile should always have been created above, for active channels
975  profileLED[fRefNum]->Fill(fHour, fAmp);
976  }
977  }
978 
979  // re-associating the branch addresses here seems to be needed for OK 'average' storage
980  fTreeLEDAvgAmpVsTime->SetBranchAddress("fRefNum", &fRefNum);
981  fTreeLEDAvgAmpVsTime->SetBranchAddress("fHour", &fHour);
982  fTreeLEDAvgAmpVsTime->SetBranchAddress("fAvgAmp", &fAvgAmp);
983  fTreeLEDAvgAmpVsTime->SetBranchAddress("fRMS", &fRMS);
984 
985  //3: fill avg tree by looping over the profiles
986  for (fRefNum = 0; fRefNum<(fgkMaxRefs*2); fRefNum++)
987  {
988  if (profileLED[fRefNum])
989  { // profile was created
990  if (profileLED[fRefNum]->GetEntries() > 0)
991  { // profile had some entries
992  for(int it=0; it<numProfBins; it++) {
993  if (profileLED[fRefNum]->GetBinEntries(it+1) > 0)
994  {
995  fAvgAmp = profileLED[fRefNum]->GetBinContent(it+1);
996  fHour = profileLED[fRefNum]->GetBinCenter(it+1);
997  fRMS = profileLED[fRefNum]->GetBinError(it+1);
998  fTreeLEDAvgAmpVsTime->Fill();
999  } // some entries for this bin
1000  } // loop over bins
1001  } // some entries for this profile
1002  } // profile exists
1003  } // loop over all possible channels
1004 
1005  // OK, we're done..
1006 
1007  return kTRUE;
1008 }
1009 
1012 //_____________________________________________________________________
1013 Bool_t AliCaloCalibSignal::DecodeChannelNum(const int chanId,
1014  int *imod, int *icol, int *irow, int *igain) const
1015 {
1016  *igain = chanId/(fModules*fColumns*fRows);
1017  *imod = (chanId/(fColumns*fRows)) % fModules;
1018  *icol = (chanId/fRows) % fColumns;
1019  *irow = chanId % fRows;
1020  return kTRUE;
1021 }
1022 
1025 //_____________________________________________________________________
1026 Bool_t AliCaloCalibSignal::DecodeRefNum(const int refId,
1027  int *imod, int *istripMod, int *igain) const
1028 {
1029  *igain = refId/(fModules*fLEDRefs);
1030  *imod = (refId/(fLEDRefs)) % fModules;
1031  *istripMod = refId % fLEDRefs;
1032  return kTRUE;
1033 }
1034 
kDetType fDetType
The detector type for this object.
int fNAcceptedEvents
Number of events accepted.
Bool_t ProcessEvent(AliRawReader *rawReader)
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
Bool_t CheckLEDRefAboveAmp(const int *AmpVal, int resultArray[]) const
static const int fgkMaxTowers
AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALCols * AliEMCALGeoParams::fgkEMCALRow...
int fStartTime
Time of first event.
int fNHighGain[fgkMaxTowers]
Number of Amp. vs. Time readings per tower.
double GetAmpCut() const
void ResetInfo()
Reset trees and counters.
bool GetUseAverage() const
void CreateTrees()
Create/initialize/setup the TTrees.
TTree * GetTreeLEDAmpVsTime() const
static Int_t GetFirstSTUDDL()
Definition: AliDAQ.h:74
kDetType
The detector types (add common enum as in AliCalibPedestal?)
int GetChannelNum(const int imod, const int icol, const int irow, const int igain) const
static const int fgkEMCALRows
Number of rows per module for EMCAL.
static const int fgkPhosRows
Number of rows per module for PHOS.
int fDownscale
To select 1 out every N (fDownscale) events.
bool fReqFractionAboveAmp
Flag to select if we should do some event selection based on amplitudes.
int fNEvents
Number of events processed.
AliCaloRawAnalyzer * fRawAnalyzer
! e.g. for sample selection for fits
AliCaloAltroMapping ** fMapping
! Altro Mapping object
static const int fgkEMCALLEDRefs
Number of LEDs (reference/monitors) per module for EMCAL; one per StripModule.
double GetHour() const
double fLatestHour
Largest fraction of hour since beginning of run, for amp vs. time graphs.
static const int fgkPhosCols
Number of columns per module for PHOS.
class for signal monitoring and calibration tools
Bool_t DecodeChannelNum(const int chanId, int *imod, int *icol, int *irow, int *igain) const
Output with the module, column, row, and gain for a given channel number.
bool GetReqFractionAboveAmp() const
void SetFittingAlgorithm(Int_t val)
int fRows
The number of rows per module.
int fRunNumber
The run number. Needs to be set by the user.
void SetParametersFromFile(const char *parameterFile)
double fReqFractionAboveAmpCutVal
Required fraction that should be above cut.
void DeleteTrees()
Delete old objects and set pointers, what was created in the constructor (TTrees) ...
TString fileName(const char *dir, int runNumber, const char *da, int i, const char *type)
Container class to hold results from fitting.
int fLEDRefs
The number of LED references/monitors per module.
Container of ALTRO information.
TTree * fTreeAvgAmpVsTime
Store channel, gain, amp, and time info, for averages.
Bool_t DecodeRefNum(const int refId, int *imod, int *istripMod, int *igain) const
Output with the module, stripModule, and gain for a given reference number.
double GetReqFractionAboveAmpCutVal() const
TTree * fTreeAmpVsTime
Store channel, gain, amp, and time info.
TString fCaloString
ID for which detector type we have.
double fAmpCut
Amplitude cut value.
static AliCaloRawAnalyzer * CreateAnalyzer(const int algo)
double fHour
Fraction of hour since beginning of run, for amp vs. time graphs, for current event.
int GetRefNum(const int imod, const int istripMod, const int igain) const
int fNRef[fgkMaxRefs *2]
Number of Amp. vs. Time readings per tower, for LED refs; *2 for both gains.
void Zero()
Set all counters to 0; not cuts etc. though.
AliCaloCalibSignal(kDetType detectorType=kEmCal)
virtual AliCaloFitResults Evaluate(const std::vector< AliCaloBunchInfo > &, UInt_t, UInt_t)=0
TTree * fTreeLEDAmpVsTime
Store channel, gain, amp, and time info, for LED reference.
bool GetReqLEDRefAboveAmpCutVal() const
void WriteParametersToFile(const char *parameterFile)
Write parameters to file.
static int fRefNum
for LED
int GetSecInAverage() const
int GetTowerNum(const int imod, const int icol, const int irow) const
int fNLowGain[fgkMaxTowers]
Number of Amp. vs. Time readings per tower, for low gain.
TTree * GetTreeAvgAmpVsTime() const
AliCaloAltroMapping * fMapping[4]
Bool_t AddInfo(const AliCaloCalibSignal *sig)
static const int fgkEMCALCols
Number of columns per module for EMCAL.
int fColumns
The number of columns per module.
double fAmpCutLEDRef
Amplitude cut value for LED reference.
int fSecInAverage
Time interval for the graph averaging.
double GetLatestHour() const
static const int fgkPhosModules
Number of modules for PHOS.
Float_t GetAmp() const
static const int fgkEMCALModules
Number of modules, 12 for EMCal + 8 for DCAL.
static int fChannelNum
for regular towers
static const int fgkNumSecInHr
Number of seconds in an hour, for the fractional hour conversion on the time graph.
TTree * GetTreeAmpVsTime() const
virtual ~AliCaloCalibSignal()
Destructor.
Bool_t CheckFractionAboveAmp(const int *AmpVal, int resultArray[]) const
bool fReqLEDRefAboveAmpCutVal
Flag to select if we should require that signal is also seen in LED Reference/Monitoring channel...
static const int fgkPhosLEDRefs
No LED monitor channels for PHOS.
int GetNAcceptedEvents() const
void res(Char_t i)
Definition: Resolution.C:2
static double fRMS
RMS.
TTree * GetTreeLEDAvgAmpVsTime() const
static double fAvgAmp
average amplitude
Bool_t Save(TString fileName)
int fModules
The number of modules.
Int_t fFittingAlgorithm
Select the fitting algorithm.
void SetIsZeroSuppressed(bool iszs=true)
TTree * fTreeLEDAvgAmpVsTime
Store channel, gain, amp, and time info, for LED reference - averages.
double GetAmpCutLEDRef() const
bool fUseAverage
Flag to average graph points into over a time interval.
static const int fgkMaxRefs
AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALLEDRefs;.
static double fAmp
amplitude