]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliCaloCalibSignal.cxx
Updated version.
[u/mrichter/AliRoot.git] / EMCAL / AliCaloCalibSignal.cxx
index e83ee0c4d775aa6a81c7ba1577f507c23119f08a..51d867073a0b5f5219c9569101d4bbeeacb2ba4a 100644 (file)
 //  fSignals->GetXXX..()
 // etc.
 //________________________________________________________________________
+#include <string>
+#include <sstream>
+#include <fstream>
 
 #include "TProfile.h"
 #include "TFile.h"
 
 #include "AliRawReader.h"
-#include "AliRawEventHeaderBase.h"
 #include "AliCaloRawStreamV3.h"
 
 //The include file
@@ -65,9 +67,11 @@ AliCaloCalibSignal::AliCaloCalibSignal(kDetType detectorType) :
   fMapping(NULL),
   fRunNumber(-1),
   fStartTime(0),
-  fAmpCut(50),
-  fReqFractionAboveAmpCutVal(0.8),
+  fAmpCut(40), // min. 40 ADC counts as default
+  fReqFractionAboveAmpCutVal(0.6), // 60% in a strip, per default
   fReqFractionAboveAmp(kTRUE),
+  fAmpCutLEDRef(100), // min. 100 ADC counts as default
+  fReqLEDRefAboveAmpCutVal(kTRUE),
   fHour(0),
   fLatestHour(0),
   fUseAverage(kTRUE),
@@ -143,6 +147,8 @@ AliCaloCalibSignal::AliCaloCalibSignal(const AliCaloCalibSignal &sig) :
   fAmpCut(sig.GetAmpCut()),
   fReqFractionAboveAmpCutVal(sig.GetReqFractionAboveAmpCutVal()),
   fReqFractionAboveAmp(sig.GetReqFractionAboveAmp()),
+  fAmpCutLEDRef(sig.GetAmpCutLEDRef()),
+  fReqLEDRefAboveAmpCutVal(sig.GetReqLEDRefAboveAmpCutVal()),
   fHour(sig.GetHour()),
   fLatestHour(sig.GetLatestHour()),
   fUseAverage(sig.GetUseAverage()),
@@ -205,7 +211,7 @@ void AliCaloCalibSignal::CreateTrees()
 
 //_____________________________________________________________________
 void AliCaloCalibSignal::ResetInfo()
-{
+{ // reset trees and counters
   Zero(); // set all counters to 0
   DeleteTrees(); // delete previous stuff
   CreateTrees(); // and create some new ones
@@ -221,7 +227,7 @@ void AliCaloCalibSignal::Zero()
   fNEvents = 0;
   fNAcceptedEvents = 0;
 
-  // Set the number of points for each tower: Amp vs. Time 
+  // Set the number of points for each tower: Amp vs. Time
   memset(fNHighGain, 0, sizeof(fNHighGain));
   memset(fNLowGain, 0, sizeof(fNLowGain));
   // and LED reference
@@ -231,28 +237,167 @@ void AliCaloCalibSignal::Zero()
 }
 
 //_____________________________________________________________________
-Bool_t AliCaloCalibSignal::CheckFractionAboveAmp(int *AmpVal, int nTotChan)
-{
-  int nAbove = 0;
+Bool_t AliCaloCalibSignal::CheckFractionAboveAmp(const int *iAmpVal, 
+                                                int resultArray[])
+{ // check fraction of towers, per column, that are above amplitude cut
+  Bool_t returnCode = false;
     
-  int TowerNum = 0;
+  int iTowerNum = 0;
+  double fraction = 0;
   for (int i = 0; i<fModules; i++) {
     for (int j = 0; j<fColumns; j++) {
+      int nAbove = 0;
       for (int k = 0; k<fRows; k++) {
-       TowerNum = GetTowerNum(i,j,k);
-       if (AmpVal[TowerNum] > fAmpCut) { 
+       iTowerNum = GetTowerNum(i,j,k);
+       if (iAmpVal[iTowerNum] > fAmpCut) { 
          nAbove++;
        }
       }
+      resultArray[i*fColumns +j] = 0; // init. to denied
+      if (nAbove > 0) {
+       fraction = (1.0*nAbove) / fRows;
+       /*
+       printf("DS mod %d col %d nAbove %d fraction %3.2f\n",
+              i, j, nAbove, fraction);
+       */
+       if (fraction > fReqFractionAboveAmpCutVal) {
+         resultArray[i*fColumns + j] = nAbove;
+         returnCode = true;
+       }  
+      }
     }
-  }
+  } // modules loop
   
-  double fraction = (1.0*nAbove) / nTotChan;
+  return returnCode;
+}
+
+
+//_____________________________________________________________________
+Bool_t AliCaloCalibSignal::CheckLEDRefAboveAmp(const int *iAmpVal, 
+                                              int resultArray[])
+{ // check which LEDRef/Mon strips are above amplitude cut
+  Bool_t returnCode = false;
+    
+  int iRefNum = 0;
+  int gain = 1; // look at high gain; this should be rather saturated usually..
+  for (int i = 0; i<fModules; i++) {
+    for (int j = 0; j<fLEDRefs; j++) {
+      iRefNum = GetRefNum(i, j, gain);
+      if (iAmpVal[iRefNum] > fAmpCutLEDRef) { 
+       resultArray[i*fLEDRefs +j] = 1; // enough signal
+       returnCode = true;
+      }
+      else {
+       resultArray[i*fLEDRefs +j] = 0; // not enough signal
+      }
+      
+      /*
+      printf("DS mod %d LEDRef %d ampVal %d\n",
+            i, j, iAmpVal[iRefNum]);
+      */
+    } // LEDRefs
+  } // modules loop
   
-  if (fraction > fReqFractionAboveAmpCutVal) {  
-    return true;
+  return returnCode;
+}
+
+// Parameter/cut handling
+//_____________________________________________________________________
+void AliCaloCalibSignal::SetParametersFromFile(const char *parameterFile)
+{ // set parameters from file
+  static const string delimitor("::");
+       
+  // open, check input file
+  ifstream in( parameterFile );
+  if( !in ) {
+    printf("in AliCaloCalibSignal::SetParametersFromFile - Using default/run_time parameters.\n");
+    return;
+  } 
+
+  // Note: this method is a bit more complicated than it really has to be
+  // - allowing for multiple entries per line, arbitrary order of the
+  // different variables etc. But I wanted to try and do this in as
+  // correct a C++ way as I could (as an exercise).
+
+  // read in
+  char readline[1024];
+  while ((in.rdstate() & ios::failbit) == 0 ) {
+    
+    // Read into the raw char array and then construct a string
+    // to do the searching
+    in.getline(readline, 1024);
+    istringstream s(readline);         
+               
+    while ( ( s.rdstate() & ios::failbit ) == 0 ) {
+                       
+      string keyValue; 
+      s >> keyValue;
+      
+      // check stream status
+      if( s.rdstate() & ios::failbit ) break;
+                       
+      // skip rest of line if comments found
+      if( keyValue.substr( 0, 2 ) == "//" ) break;
+                       
+      // look for "::" in keyValue pair
+      size_t position = keyValue.find( delimitor );
+      if( position == string::npos ) {
+       printf("wrong format for key::value pair: %s\n", keyValue.c_str());
+      }
+                               
+      // split keyValue pair
+      string key( keyValue.substr( 0, position ) );
+      string value( keyValue.substr( position+delimitor.size(), 
+                                     keyValue.size()-delimitor.size() ) );
+                       
+      // check value does not contain a new delimitor
+      if( value.find( delimitor ) != string::npos ) {
+       printf("wrong format for key::value pair: %s\n", keyValue.c_str());
+      }
+      
+      // debug: check key value pair
+      // printf("AliCaloCalibSignal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
+
+      // if the key matches with something we expect, we assign the new value
+      if ( (key == "fAmpCut") || (key == "fReqFractionAboveAmpCutVal") ||
+          (key == "fAmpCutLEDRef") || (key == "fSecInAverage") ) {
+       istringstream iss(value);
+       printf("AliCaloCalibSignal::SetParametersFromFile - key %s value %s\n", key.c_str(), value.c_str());
+
+       if (key == "fAmpCut") { 
+         iss >> fAmpCut; 
+       }
+       else if (key == "fReqFractionAboveAmpCutVal") { 
+         iss >> fReqFractionAboveAmpCutVal; 
+       }
+       else if (key == "fAmpCutLEDRef") { 
+         iss >> fAmpCutLEDRef; 
+       }
+       else if (key == "fSecInAverage") { 
+         iss >> fSecInAverage; 
+       }
+      } // some match found/expected
+
+    }          
   }
-  else return false;
+
+  in.close();
+  return;      
+}
+
+//_____________________________________________________________________
+void AliCaloCalibSignal::WriteParametersToFile(const char *parameterFile)
+{ // write parameters to file
+  static const string delimitor("::");
+  ofstream out( parameterFile );
+  out << "// " << parameterFile << endl;
+  out << "fAmpCut" << "::" << fAmpCut << endl;
+  out << "fReqFractionAboveAmpCutVal" << "::" << fReqFractionAboveAmpCutVal << endl;
+  out << "fAmpCutLEDRef" << "::" << fAmpCutLEDRef << endl;
+  out << "fSecInAverage" << "::" << fSecInAverage << endl;
+
+  out.close();
+  return;
 }
 
 //_____________________________________________________________________
@@ -324,13 +469,15 @@ Bool_t AliCaloCalibSignal::AddInfo(const AliCaloCalibSignal *sig)
 Bool_t AliCaloCalibSignal::ProcessEvent(AliRawReader *rawReader)
 {
   // if fMapping is NULL the rawstream will crate its own mapping
-  AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);
-
-  return ProcessEvent( &rawStream, (AliRawEventHeaderBase*)rawReader->GetEventHeader() );
+  AliCaloRawStreamV3 rawStream(rawReader, fCaloString, (AliAltroMapping**)fMapping);  
+  if (fDetType == kEmCal) {
+    rawReader->Select("EMCAL", 0, AliEMCALGeoParams::fgkLastAltroDDL) ; //select EMCAL DDL range 
+  }
+  return ProcessEvent( &rawStream, rawReader->GetTimestamp() );
 }
 
 //_____________________________________________________________________
-Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStreamV3 *in, AliRawEventHeaderBase *aliHeader)
+Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStreamV3 *in, UInt_t Timestamp)
 { 
   // Method to process=analyze one event in the data stream
   if (!in) return kFALSE; //Return right away if there's a null pointer
@@ -338,24 +485,25 @@ Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStreamV3 *in, AliRawEventHeade
   fNEvents++; // one more event
 
   // use maximum numbers to set array sizes
-  int AmpValHighGain[fgkMaxTowers];
-  int AmpValLowGain[fgkMaxTowers];
-  memset(AmpValHighGain, 0, sizeof(AmpValHighGain));
-  memset(AmpValLowGain, 0, sizeof(AmpValLowGain));
+  int iAmpValHighGain[fgkMaxTowers];
+  int iAmpValLowGain[fgkMaxTowers];
+  memset(iAmpValHighGain, 0, sizeof(iAmpValHighGain));
+  memset(iAmpValLowGain, 0, sizeof(iAmpValLowGain));
 
   // also for LED reference
-  int LEDAmpVal[fgkMaxRefs * 2]; // factor 2 is for the two gain values
-  memset(LEDAmpVal, 0, sizeof(LEDAmpVal));
+  int iLEDAmpVal[fgkMaxRefs * 2]; // factor 2 is for the two gain values
+  memset(iLEDAmpVal, 0, sizeof(iLEDAmpVal));
 
   int sample; // temporary value
   int gain = 0; // high or low gain
   
-  // Number of Low and High gain channels for this event:
+  // Number of Low and High gain, and LED Ref, channels for this event:
   int nLowChan = 0; 
   int nHighChan = 0; 
+  int nLEDRefChan = 0;
 
-  int TowerNum = 0; // array index for regular towers
-  int RefNum = 0; // array index for LED references
+  int iTowerNum = 0; // array index for regular towers
+  int iRefNum = 0; // array index for LED references
 
   // loop first to get the fraction of channels with amplitudes above cut
 
@@ -403,23 +551,24 @@ Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStreamV3 *in, AliRawEventHeade
 
       if ( in->IsHighGain() || in->IsLowGain() ) { // regular tower
        // get tower number for AmpVal array
-       TowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow()); 
+       iTowerNum = GetTowerNum(arrayPos, in->GetColumn(), in->GetRow()); 
 
        if (gain == 0) {
          // fill amplitude into the array         
-         AmpValLowGain[TowerNum] = max - min;
+         iAmpValLowGain[iTowerNum] = max - min;
          nLowChan++;
        } 
        else if (gain==1) {//fill the high gain ones
          // fill amplitude into the array
-         AmpValHighGain[TowerNum] = max - min;
+         iAmpValHighGain[iTowerNum] = max - min;
          nHighChan++;
        }//end if gain
       } // regular tower
       else if ( in->IsLEDMonData() ) { // LED ref.; 
        // strip # is coded is 'column' in the channel maps 
-       RefNum = GetRefNum(arrayPos, in->GetColumn(), gain); 
-       LEDAmpVal[RefNum] = max - min;
+       iRefNum = GetRefNum(arrayPos, in->GetColumn(), gain); 
+       iLEDAmpVal[iRefNum] = max - min;
+       nLEDRefChan++;
       } // end of LED ref
 
       } // nsamples>0 check, some data found for this channel; not only trailer/header      
@@ -429,11 +578,30 @@ Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStreamV3 *in, AliRawEventHeade
   
   in->Reset(); // just in case the next customer forgets to check if the stream was reset..
 
-  // now check if it was a led event, only use high gain (that should be sufficient)
+  // now check if it was an LED event, using the LED Reference info per strip
+
+  // by default all columns are accepted (init check to > 0)
+  int checkResultArray[AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALCols];
+  for (int ia=0; ia<(AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALCols); ia++) { 
+    checkResultArray[ia] = 1; 
+  }
   if (fReqFractionAboveAmp) {
     bool ok = false;
     if (nHighChan > 0) { 
-      ok = CheckFractionAboveAmp(AmpValHighGain, nHighChan); 
+      ok = CheckFractionAboveAmp(iAmpValHighGain, checkResultArray); 
+    }
+    if (!ok) return false; // skip event
+  }
+
+  // by default all columns are accepted (init check to > 0)
+  int checkResultArrayLEDRef[AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALLEDRefs];
+  for (int ia=0; ia<(AliEMCALGeoParams::fgkEMCALModules * AliEMCALGeoParams::fgkEMCALLEDRefs); ia++) { 
+    checkResultArrayLEDRef[ia] = 1; 
+  }
+  if (fReqLEDRefAboveAmpCutVal) {
+    bool ok = false;
+    if (nLEDRefChan > 0) { 
+      ok = CheckLEDRefAboveAmp(iLEDAmpVal, checkResultArrayLEDRef); 
     }
     if (!ok) return false; // skip event
   }
@@ -441,46 +609,60 @@ Bool_t AliCaloCalibSignal::ProcessEvent(AliCaloRawStreamV3 *in, AliRawEventHeade
   fNAcceptedEvents++; // one more event accepted
 
   if (fStartTime == 0) { // if start-timestamp wasn't set,we'll pick it up from the first event we encounter
-    fStartTime = aliHeader->Get("Timestamp");
+    fStartTime = Timestamp;
   }
 
-  fHour = (aliHeader->Get("Timestamp")-fStartTime)/(double)fgkNumSecInHr;
+  fHour = (Timestamp - fStartTime)/(double)fgkNumSecInHr;
   if (fLatestHour < fHour) {
     fLatestHour = fHour; 
   }
   
   // it is a led event, now fill TTree
+  // We also do the activity check for LEDRefs/Strips, but need to translate between column
+  // and strip indices for that; based on these relations: 
+  // iStrip = AliEMCALGeoParams::GetStripModule(iSM, iCol);
+  // iStrip = (iSM%2==0) ? iCol/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iCol/2;
+  // which leads to
+  // iColFirst = (iSM%2==0) ? iStrip*2 : (AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - iStrip)*2;
+
   for(int i=0; i<fModules; i++){
-    for(int j=0; j<fColumns; j++){
+    for(int j=0; j<fColumns; j++) {
+      int iStrip = (i%2==0) ? j/2 : AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - j/2;
+      if (checkResultArray[i*fColumns + j]>0  && checkResultArrayLEDRef[i*fLEDRefs + iStrip]>0) { // column passed check 
       for(int k=0; k<fRows; k++){
        
-       TowerNum = GetTowerNum(i, j, k); 
+       iTowerNum = GetTowerNum(i, j, k); 
 
-       if(AmpValHighGain[TowerNum]) {
-         fAmp = AmpValHighGain[TowerNum];
+       if(iAmpValHighGain[iTowerNum]) {
+         fAmp = iAmpValHighGain[iTowerNum];
          fChannelNum = GetChannelNum(i,j,k,1);
-         fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValHighGain[TowerNum]);
-         fNHighGain[TowerNum]++;
+         fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValHighGain[iTowerNum]);
+         fNHighGain[iTowerNum]++;
        }
-       if(AmpValLowGain[TowerNum]) {
-         fAmp = AmpValLowGain[TowerNum];
+       if(iAmpValLowGain[iTowerNum]) {
+         fAmp = iAmpValLowGain[iTowerNum];
          fChannelNum = GetChannelNum(i,j,k,0);
-         fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValLowGain[TowerNum]);
-         fNLowGain[TowerNum]++;
+         fTreeAmpVsTime->Fill();//fChannelNum,fHour,AmpValLowGain[iTowerNum]);
+         fNLowGain[iTowerNum]++;
        }
       } // rows
+      } // column passed check, and LED Ref for strip passed check (if any)
     } // columns
 
     // also LED refs
     for(int j=0; j<fLEDRefs; j++){
-      for (gain=0; gain<2; gain++) {
-       fRefNum = GetRefNum(i, j, gain); 
-       if (LEDAmpVal[fRefNum]) {
-         fAmp = LEDAmpVal[fRefNum];
-         fTreeLEDAmpVsTime->Fill();//fRefNum,fHour,fAmp);
-         fNRef[fRefNum]++;
-       }
-      }
+      int iColFirst = (i%2==0) ? j*2 : (AliEMCALGeoParams::fgkEMCALLEDRefs - 1 - j)*2; //CHECKME!!!
+      if ( ((checkResultArray[i*fColumns + iColFirst]>0) || (checkResultArray[i*fColumns + iColFirst + 1]>0)) && // at least one column in strip passed check 
+          (checkResultArrayLEDRef[i*fLEDRefs + j]>0) ) { // and LED Ref passed checks
+       for (gain=0; gain<2; gain++) {
+         fRefNum = GetRefNum(i, j, gain); 
+         if (iLEDAmpVal[fRefNum]) {
+           fAmp = iLEDAmpVal[fRefNum];
+           fTreeLEDAmpVsTime->Fill();//fRefNum,fHour,fAmp);
+           fNRef[fRefNum]++;
+         }
+       } // gain
+      } // at least one column in strip passed check, and LED Ref passed check (if any) 
     }
 
   } // modules
@@ -544,22 +726,22 @@ Bool_t AliCaloCalibSignal::Analyze()
   memset(profile, 0, sizeof(profile));
 
   char name[200]; // for profile id and title
-  int TowerNum = 0;
+  int iTowerNum = 0;
 
   for (int i = 0; i<fModules; i++) {
     for (int ic=0; ic<fColumns; ic++){
       for (int ir=0; ir<fRows; ir++) {
 
-       TowerNum = GetTowerNum(i, ic, ir);
+       iTowerNum = GetTowerNum(i, ic, ir);
        // high gain
-       if (fNHighGain[TowerNum] > 0) {
+       if (fNHighGain[iTowerNum] > 0) {
          fChannelNum = GetChannelNum(i, ic, ir, 1); 
          sprintf(name, "profileChan%d", fChannelNum);
          profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
        }
 
        // same for low gain
-       if (fNLowGain[TowerNum] > 0) {
+       if (fNLowGain[iTowerNum] > 0) {
          fChannelNum = GetChannelNum(i, ic, ir, 0); 
          sprintf(name, "profileChan%d", fChannelNum);
          profile[fChannelNum] = new TProfile(name, name, numProfBins, timeMin, timeMax, "s");
@@ -664,3 +846,24 @@ Bool_t AliCaloCalibSignal::Analyze()
 
   return kTRUE;
 }
+
+//_____________________________________________________________________
+Bool_t AliCaloCalibSignal::DecodeChannelNum(const int chanId, 
+                                           int *imod, int *icol, int *irow, int *igain) const  
+{ // return the module, column, row, and gain for a given channel number
+  *igain = chanId/(fModules*fColumns*fRows);
+  *imod = (chanId/(fColumns*fRows)) % fModules;
+  *icol = (chanId/fRows) % fColumns;
+  *irow = chanId % fRows;
+  return kTRUE;
+} 
+
+//_____________________________________________________________________
+Bool_t AliCaloCalibSignal::DecodeRefNum(const int refId, 
+                                       int *imod, int *istripMod, int *igain) const 
+{ // return the module, stripModule, and gain for a given reference number
+  *igain = refId/(fModules*fLEDRefs);
+  *imod = (refId/(fLEDRefs)) % fModules;
+  *istripMod = refId % fLEDRefs;
+  return kTRUE;
+}