]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCCalibCE.cxx
Use debug stream only if requested
[u/mrichter/AliRoot.git] / TPC / AliTPCCalibCE.cxx
index d4e425457bff06ce3bdb3dbab36b80b3eaf6c853..6e1c4c81f36e433839f43b9d07e805ae3e045653 100644 (file)
@@ -275,6 +275,7 @@ END_HTML */
 #include <TFile.h>
 
 //AliRoot includes
+#include "AliLog.h"
 #include "AliRawReader.h"
 #include "AliRawReaderRoot.h"
 #include "AliRawReaderDate.h"
@@ -308,8 +309,13 @@ AliTPCCalibCE::AliTPCCalibCE() :
     fNbinsRMS(100),
     fXminRMS(0.1),
     fXmaxRMS(5.1),
+    fPeakMinus(2),
+    fPeakPlus(3),
+    fNoiseThresholdMax(5.),
+    fNoiseThresholdSum(8.),
+    fIsZeroSuppressed(kFALSE),
     fLastSector(-1),
-    fOldRCUformat(kTRUE),
+    fSecRejectRatio(.4),
     fROC(AliTPCROC::Instance()),
     fMapping(NULL),
     fParam(new AliTPCParam),
@@ -318,12 +324,16 @@ AliTPCCalibCE::AliTPCCalibCE() :
     fPedestalROC(0x0),
     fPadNoiseROC(0x0),
     fCalRocArrayT0(72),
+    fCalRocArrayT0Err(72),
     fCalRocArrayQ(72),
     fCalRocArrayRMS(72),
     fCalRocArrayOutliers(72),
     fHistoQArray(72),
     fHistoT0Array(72),
     fHistoRMSArray(72),
+    fMeanT0rms(0),
+    fMeanQrms(0),
+    fMeanRMSrms(0),
     fHistoTmean(72),
     fParamArrayEventPol1(72),
     fParamArrayEventPol2(72),
@@ -360,6 +370,7 @@ AliTPCCalibCE::AliTPCCalibCE() :
     // AliTPCSignal default constructor
     //
 //    fHTime0 = new TH1F("hTime0Event","hTime0Event",(fLastTimeBin-fFirstTimeBin)*10,fFirstTimeBin,fLastTimeBin);
+    fParam->Update();
 }
 //_____________________________________________________________________
 AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
@@ -375,8 +386,13 @@ AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
     fNbinsRMS(sig.fNbinsRMS),
     fXminRMS(sig.fXminRMS),
     fXmaxRMS(sig.fXmaxRMS),
+    fPeakMinus(sig.fPeakMinus),
+    fPeakPlus(sig.fPeakPlus),
+    fNoiseThresholdMax(sig.fNoiseThresholdMax),
+    fNoiseThresholdSum(sig.fNoiseThresholdSum),
+    fIsZeroSuppressed(sig.fIsZeroSuppressed),
     fLastSector(-1),
-    fOldRCUformat(kTRUE),
+    fSecRejectRatio(.4),
     fROC(AliTPCROC::Instance()),
     fMapping(NULL),
     fParam(new AliTPCParam),
@@ -385,12 +401,16 @@ AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
     fPedestalROC(0x0),
     fPadNoiseROC(0x0),
     fCalRocArrayT0(72),
+    fCalRocArrayT0Err(72),
     fCalRocArrayQ(72),
     fCalRocArrayRMS(72),
     fCalRocArrayOutliers(72),
     fHistoQArray(72),
     fHistoT0Array(72),
     fHistoRMSArray(72),
+    fMeanT0rms(sig.fMeanT0rms),
+    fMeanQrms(sig.fMeanQrms),
+    fMeanRMSrms(sig.fMeanRMSrms),
     fHistoTmean(72),
     fParamArrayEventPol1(72),
     fParamArrayEventPol2(72),
@@ -503,6 +523,7 @@ AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
     fVEventTime.SetElements(sig.fVEventTime.GetMatrixArray());
     fVEventNumber.SetElements(sig.fVEventNumber.GetMatrixArray());
 
+    fParam->Update();
 }
 //_____________________________________________________________________
 AliTPCCalibCE& AliTPCCalibCE::operator = (const  AliTPCCalibCE &source)
@@ -523,6 +544,7 @@ AliTPCCalibCE::~AliTPCCalibCE()
     //
 
     fCalRocArrayT0.Delete();
+    fCalRocArrayT0Err.Delete();
     fCalRocArrayQ.Delete();
     fCalRocArrayRMS.Delete();
     fCalRocArrayOutliers.Delete();
@@ -561,6 +583,11 @@ Int_t AliTPCCalibCE::Update(const Int_t icsector,
     // assumes that it is looped over consecutive time bins of one pad
     //
 
+    //temp
+//    if (icsector<36) return 0;
+//    if (icsector%36>17) return 0;
+
+
   if (icRow<0) return 0;
   if (icPad<0) return 0;
   if (icTimeBin<0) return 0;
@@ -610,7 +637,7 @@ void AliTPCCalibCE::FindPedestal(Float_t part)
        }
 
        if ( fPedestalROC&&fPadNoiseROC ){
-           fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
+           fPadPedestal = fPedestalROC->GetValue(fCurrentChannel)*(Float_t)(!fIsZeroSuppressed);
            fPadNoise    = fPadNoiseROC->GetValue(fCurrentChannel);
             noPedestal   = kFALSE;
        }
@@ -620,6 +647,9 @@ void AliTPCCalibCE::FindPedestal(Float_t part)
     //if we are not running with pedestal database, or for the current sector there is no information
     //available, calculate the pedestal and noise on the fly
     if ( noPedestal ) {
+       fPadPedestal = 0;
+       fPadNoise    = 0;
+        if ( fIsZeroSuppressed ) return;
        const Int_t kPedMax = 100;  //maximum pedestal value
        Float_t  max    =  0;
        Float_t  maxPos =  0;
@@ -665,8 +695,6 @@ void AliTPCCalibCE::FindPedestal(Float_t part)
                rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
            }
        }
-       fPadPedestal = 0;
-       fPadNoise    = 0;
        if ( count > 0 ) {
            mean/=count;
            rms    = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
@@ -686,7 +714,7 @@ void AliTPCCalibCE::FindCESignal(TVectorD &param, Float_t &qSum, const TVectorF
 
     Float_t ceQmax  =0, ceQsum=0, ceTime=0, ceRMS=0;
     Int_t   cemaxpos       = 0;
-    Float_t ceSumThreshold = 8.*fPadNoise;  // threshold for the signal sum
+    Float_t ceSumThreshold = fNoiseThresholdSum*fPadNoise;  // threshold for the signal sum
     const Int_t    kCemin  = 4;             // range for the analysis of the ce signal +- channels from the peak
     const Int_t    kCemax  = 7;
 
@@ -761,12 +789,12 @@ void AliTPCCalibCE::FindLocalMaxima(TVectorF &maxima)
     //
     // Find local maxima on the pad signal and Histogram them
     //
-  Float_t ceThreshold = 5.*TMath::Max(fPadNoise,Float_t(1.));  // threshold for the signal
+  Float_t ceThreshold = fNoiseThresholdMax*TMath::Max(fPadNoise,Float_t(1.));  // threshold for the signal
     Int_t   count       = 0;
-    Int_t   tminus      = 2;
-    Int_t   tplus       = 3;
-    for (Int_t i=fLastTimeBin-tplus-1; i>=fFirstTimeBin+tminus; --i){
-       if ( (fPadSignal[i]-fPadPedestal)>ceThreshold && IsPeak(i,tminus,tplus) ){
+//    Int_t   tminus      = 2;
+//    Int_t   tplus       = 3;
+    for (Int_t i=fLastTimeBin-fPeakPlus-1; i>=fFirstTimeBin+fPeakMinus; --i){
+       if ( (fPadSignal[i]-fPadPedestal)>ceThreshold && IsPeak(i,fPeakMinus,fPeakPlus) ){
          if (count<maxima.GetNrows()){
            maxima.GetMatrixArray()[count++]=i;
            GetHistoTmean(fCurrentSector,kTRUE)->Fill(i);
@@ -788,6 +816,8 @@ void AliTPCCalibCE::ProcessPad()
     FindLocalMaxima(maxima);
     if ( (fNevents == 0) || (fOldRunNumber!=fRunNumber) ) return;  // return because we don't have Time0 info for the CE yet
 
+    if ( !GetTMeanEvents(fCurrentSector) ) return; //return if we don't have time 0 info, eg if only one side has laser
+
     TVectorD param(3);
     Float_t  qSum;
     FindCESignal(param, qSum, maxima);
@@ -828,6 +858,8 @@ void AliTPCCalibCE::EndEvent()
     //check if last pad has allready been processed, if not do so
     if ( fMaxTimeBin>-1 ) ProcessPad();
 
+    AliDebug(3, Form("EndEvent() - Start; Event: %05d", fNevents));
+
     TVectorD param(3);
     TMatrixD dummy(3,3);
 //    TVectorF vMeanTime(72);
@@ -849,11 +881,19 @@ void AliTPCCalibCE::EndEvent()
        time0Side[1]/=time0SideCount[1];
     // end find time0 offset
 
+    Int_t nSecMeanT=0;
     //loop over all ROCs, fill CE Time histogram corrected for the mean Time0 of each ROC
     for ( Int_t iSec = 0; iSec<72; ++iSec ){
+       AliDebug(4,Form("Processing sector '%02d'\n",iSec));
       //find median and then calculate the mean around it
        TH1S *hMeanT    = GetHistoTmean(iSec); //histogram with local maxima position information
        if ( !hMeanT ) continue;
+        //continue if not enough data is filled in the meanT histogram. This is the case if we do not have a laser event.
+       if ( hMeanT->GetEntries() < fROC->GetNChannels(iSec)*fSecRejectRatio ){
+           hMeanT->Reset();
+           AliDebug(3,Form("Skipping sec. '%02d': Not enough statistics\n",iSec));
+          continue;
+        }
 
        Double_t entries = hMeanT->GetEntries();
        Double_t sum     = 0;
@@ -861,7 +901,7 @@ void AliTPCCalibCE::EndEvent()
         Int_t ibin=0;
        for ( ibin=0; ibin<hMeanT->GetNbinsX(); ++ibin){
            sum+=arr[ibin];
-            if ( sum>=(entries/2) ) break;
+            if ( sum>=(entries/2.) ) break;
        }
        Int_t delta = 4;
         Int_t firstBin = fFirstTimeBin+ibin-delta;
@@ -877,8 +917,9 @@ void AliTPCCalibCE::EndEvent()
             vMeanTime->ResizeTo(vSize+100);
 
        vMeanTime->GetMatrixArray()[fNevents]=median;
+        nSecMeanT++;
       // end find median
-
+        
        TVectorF *vTimes = GetPadTimesEvent(iSec);
        if ( !vTimes ) continue;                     //continue if no time information for this sector is available
 
@@ -972,7 +1013,6 @@ void AliTPCCalibCE::EndEvent()
            }
            //-----------------------------  Debug end  ------------------------------
        }// end channel loop
-       hMeanT->Reset();
 
        TVectorD paramPol1(3);
        TVectorD paramPol2(6);
@@ -993,7 +1033,6 @@ void AliTPCCalibCE::EndEvent()
            GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
            GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
        }
-//     printf("events: %d -- size: %d\n",fNevents,GetParamArrayPol1(iSec)->GetSize());
 
        //-------------------------------  Debug start  ------------------------------
        if ( fDebugLevel>0 ){
@@ -1019,8 +1058,12 @@ void AliTPCCalibCE::EndEvent()
                "\n";
        }
        //-------------------------------  Debug end  ------------------------------
+        hMeanT->Reset();
     }// end sector loop
-
+    //return if no sector has a valid mean time
+    if ( nSecMeanT == 0 ) return;
+    
+    
 //    fTMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanTime),fNevents);
 //    fQMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanQ),fNevents);
     if ( fVEventTime.GetNrows() < fNevents+1 ) {
@@ -1035,6 +1078,7 @@ void AliTPCCalibCE::EndEvent()
 
     delete calIroc;
     delete calOroc;
+    AliDebug(3, Form("EndEvent() - End; Event: %05d", fNevents));
 }
 //_____________________________________________________________________
 Bool_t AliTPCCalibCE::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
@@ -1045,24 +1089,24 @@ Bool_t AliTPCCalibCE::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
   ResetEvent();
   Bool_t withInput = kFALSE;
   while ( rawStreamFast->NextDDL() ){
-      while ( rawStreamFast->NextChannel() ){
-         Int_t isector  = rawStreamFast->GetSector();                       //  current sector
-         Int_t iRow     = rawStreamFast->GetRow();                          //  current row
-         Int_t iPad     = rawStreamFast->GetPad();                          //  current pad
-         Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
-          Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
-
-         while ( rawStreamFast->NextBunch() ){
-             for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
-                 Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin];
-                 Update(isector,iRow,iPad,iTimeBin+1,signal);
-                 withInput = kTRUE;
-             }
-         }
+    while ( rawStreamFast->NextChannel() ){
+      Int_t isector  = rawStreamFast->GetSector();                       //  current sector
+      Int_t iRow     = rawStreamFast->GetRow();                          //  current row
+      Int_t iPad     = rawStreamFast->GetPad();                          //  current pad
+
+      while ( rawStreamFast->NextBunch() ){
+        Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
+        Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
+        for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
+          Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin];
+         Update(isector,iRow,iPad,iTimeBin+1,signal);
+         withInput = kTRUE;
+       }
       }
+    }
   }
   if (withInput){
-      EndEvent();
+    EndEvent();
   }
   return withInput;
 }
@@ -1096,8 +1140,6 @@ Bool_t AliTPCCalibCE::ProcessEvent(AliTPCRawStream *rawStream)
   // The Function 'SetTimeStamp' should be called for each event to set the event time stamp!!!
   //
 
-  rawStream->SetOldRCUFormat(fOldRCUformat);
-
   ResetEvent();
 
   Bool_t withInput = kFALSE;
@@ -1226,7 +1268,7 @@ TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
     sprintf(name,"hCalib%s%.2d",type,sector);
     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
 
-    // new histogram with calib information. One value for each pad!
+    // new histogram with calib information. One value for each pad!
     TH1S* hist = new TH1S(name,title,
                          fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
     hist->SetDirectory(0);
@@ -1342,13 +1384,23 @@ AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t forc
 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force)
 {
     //
-    // return pointer to Carge ROC Calibration
+    // return pointer to Time 0 ROC Calibration
     // if force is true create a new histogram if it doesn't exist allready
     //
     TObjArray *arr = &fCalRocArrayT0;
     return GetCalRoc(sector, arr, force);
 }
 //_____________________________________________________________________
+AliTPCCalROC* AliTPCCalibCE::GetCalRocT0Err(Int_t sector, Bool_t force)
+{
+    //
+    // return pointer to the error of Time 0 ROC Calibration
+    // if force is true create a new histogram if it doesn't exist allready
+    //
+    TObjArray *arr = &fCalRocArrayT0Err;
+    return GetCalRoc(sector, arr, force);
+}
+//_____________________________________________________________________
 AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force)
 {
     //
@@ -1652,14 +1704,20 @@ void AliTPCCalibCE::Analyse()
     TVectorD paramRMS(3);
     TMatrixD dummy(3,3);
 
+    Float_t channelCounter=0;
+    fMeanT0rms=0;
+    fMeanQrms=0;
+    fMeanRMSrms=0;
+
     for (Int_t iSec=0; iSec<72; ++iSec){
        TH2S *hT0 = GetHistoT0(iSec);
         if (!hT0 ) continue;
 
-       AliTPCCalROC *rocQ   = GetCalRocQ  (iSec,kTRUE);
-       AliTPCCalROC *rocT0  = GetCalRocT0 (iSec,kTRUE);
-       AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
-        AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
+       AliTPCCalROC *rocQ     = GetCalRocQ  (iSec,kTRUE);
+       AliTPCCalROC *rocT0    = GetCalRocT0 (iSec,kTRUE);
+       AliTPCCalROC *rocT0Err = GetCalRocT0Err (iSec,kTRUE);
+       AliTPCCalROC *rocRMS   = GetCalRocRMS(iSec,kTRUE);
+        AliTPCCalROC *rocOut   = GetCalRocOutliers(iSec,kTRUE);
 
        TH2S *hQ   = GetHistoQ(iSec);
        TH2S *hRMS = GetHistoRMS(iSec);
@@ -1682,18 +1740,22 @@ void AliTPCCalibCE::Analyse()
            Float_t cogTime0 = -1000;
            Float_t cogQ     = -1000;
            Float_t cogRMS   = -1000;
-            Float_t cogOut   = 0;
+           Float_t cogOut   = 0;
+            Float_t rms      = 0;
+            Float_t rmsT0    = 0;
 
 
            Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
            Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
            Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
 
-           cogQ     = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
-           cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
-            cogRMS   = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
-
-
+           cogQ     = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&rms);
+            fMeanQrms+=rms;
+           cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&rmsT0);
+            fMeanT0rms+=rmsT0;
+            cogRMS   = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&rms);
+            fMeanRMSrms+=rms;
+            channelCounter++;
 
            /*
              //outlier specifications
@@ -1706,6 +1768,7 @@ void AliTPCCalibCE::Analyse()
 */
                    rocQ->SetValue(iChannel, cogQ*cogQ);
            rocT0->SetValue(iChannel, cogTime0);
+           rocT0Err->SetValue(iChannel, rmsT0);
            rocRMS->SetValue(iChannel, cogRMS);
            rocOut->SetValue(iChannel, cogOut);
 
@@ -1739,6 +1802,11 @@ void AliTPCCalibCE::Analyse()
        }
 
     }
+    if ( channelCounter>0 ){
+       fMeanT0rms/=channelCounter;
+       fMeanQrms/=channelCounter;
+       fMeanRMSrms/=channelCounter;
+    }
     if ( fDebugStreamer ) fDebugStreamer->GetFile()->Write();
 //    delete fDebugStreamer;
 //    fDebugStreamer = 0x0;