]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDclusterizer.cxx
streamed selection of the backup track for refit and moved to
[u/mrichter/AliRoot.git] / TRD / AliTRDclusterizer.cxx
index 380596e483704c7eb33268aa3913ca5aed7b2575..191ec02b24d2e09b11c04d812b49b9f30319b0f8 100644 (file)
@@ -85,6 +85,7 @@ AliTRDclusterizer::AliTRDclusterizer(const AliTRDReconstructor *const rec)
   ,firstClusterROC(0)
   ,fNoOfClusters(0)
   ,fBaseline(0)
+  ,fRawStream(NULL)
 {
   //
   // AliTRDclusterizer default constructor
@@ -144,6 +145,7 @@ AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title, co
   ,firstClusterROC(0)
   ,fNoOfClusters(0)
   ,fBaseline(0)
+  ,fRawStream(NULL)
 {
   //
   // AliTRDclusterizer constructor
@@ -198,6 +200,7 @@ AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c)
   ,firstClusterROC(0)
   ,fNoOfClusters(0)
   ,fBaseline(0)
+  ,fRawStream(NULL)
 {
   //
   // AliTRDclusterizer copy constructor
@@ -241,7 +244,12 @@ AliTRDclusterizer::~AliTRDclusterizer()
 
   if (fTransform){
     delete fTransform;
-    fTransform     = NULL;
+    fTransform = NULL;
+  }
+
+  if (fRawStream){
+    delete fRawStream;
+    fRawStream = NULL;
   }
 
 }
@@ -296,6 +304,7 @@ void AliTRDclusterizer::Copy(TObject &c) const
   ((AliTRDclusterizer &) c).firstClusterROC= 0;
   ((AliTRDclusterizer &) c).fNoOfClusters  = 0;
   ((AliTRDclusterizer &) c).fBaseline      = 0;
+  ((AliTRDclusterizer &) c).fRawStream     = NULL;
   
 }
 
@@ -444,23 +453,29 @@ Bool_t AliTRDclusterizer::WriteClusters(Int_t det)
       ioArray->AddLast(c);
     }
     fClusterTree->Fill();
+    ioArray->Clear();
   } else {
-    
-    Int_t detOld = -1;
+    Int_t detOld = -1, nw(0);
     for (Int_t i = 0; i < nRecPoints; i++) {
       AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
       if(c->GetDetector() != detOld){
+        nw += ioArray->GetEntriesFast();
         fClusterTree->Fill();
         ioArray->Clear();
         detOld = c->GetDetector();
       } 
       ioArray->AddLast(c);
     }
+    if(ioArray->GetEntriesFast()){
+      nw += ioArray->GetEntriesFast();
+      fClusterTree->Fill();
+      ioArray->Clear();
+    }
+    AliDebug(2, Form("Clusters FOUND[%d] WRITTEN[%d] STATUS[%s]", nRecPoints, nw, nw==nRecPoints?"OK":"FAILED"));
   }
   delete ioArray;
 
   return kTRUE;  
-
 }
 
 //_____________________________________________________________________________
@@ -596,10 +611,8 @@ Bool_t AliTRDclusterizer::MakeClusters()
     //  ResetRecPoints();
     //}
         
-    // No compress just remove
-    fDigitsManager->RemoveDigits(i);
-    fDigitsManager->RemoveDictionaries(i);      
-    fDigitsManager->ClearIndexes(i);  
+    // Clear arrays of this chamber, to prepare for next event
+    fDigitsManager->ClearArrays(i);
   }
   
   if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
@@ -646,18 +659,24 @@ Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
     fTrackletContainer[1] = new UInt_t[kTrackletChmb]; 
   }
 
-  AliTRDrawStreamBase *input = AliTRDrawStreamBase::GetRawStream(rawReader);
-  if(fReconstructor->IsHLT())
-    input->SetSharedPadReadout(kFALSE);
+  if(!fRawStream)
+    fRawStream = AliTRDrawStreamBase::GetRawStream(rawReader);
+  else
+    fRawStream->SetReader(rawReader);
+
+  SetBit(kHLT, fReconstructor->IsHLT());
+
+  if(TestBit(kHLT)){
+    fRawStream->SetSharedPadReadout(kFALSE);
+    fRawStream->SetNoErrorWarning();
+  }
 
-  AliInfo(Form("Stream version: %s", input->IsA()->GetName()));
+  AliDebug(1,Form("Stream version: %s", fRawStream->IsA()->GetName()));
   
   Int_t det    = 0;
-  while ((det = input->NextChamber(fDigitsManager,fTrackletContainer)) >= 0){
-    Bool_t iclusterBranch = kFALSE;
-    if (fDigitsManager->GetIndexes(det)->HasEntry()){
-      iclusterBranch = MakeClusters(det);
-    }
+  while ((det = fRawStream->NextChamber(fDigitsManager,fTrackletContainer)) >= 0){
+    if (fDigitsManager->GetIndexes(det)->HasEntry())
+      MakeClusters(det);
 
     fDigitsManager->ClearArrays(det);
 
@@ -677,11 +696,10 @@ Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
   if(!TestBit(knewDM)){
     delete fDigitsManager;
     fDigitsManager = NULL;
+    delete fRawStream;
+    fRawStream = NULL;
   }
 
-  delete input;
-  input = NULL;
-
   AliInfo(Form("Number of found clusters : %d", fNoOfClusters)); 
   return kTRUE;
 
@@ -742,7 +760,7 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
 
   // Get the digits
   fDigits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod
-  fBaseline = fDigitsManager->GetDigitsParam()->GetADCbaseline();
+  fBaseline = fDigitsManager->GetDigitsParam()->GetADCbaseline(det);
   
   // This is to take care of switched off super modules
   if (!fDigits->HasData()) return kFALSE;
@@ -753,7 +771,7 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
     return kFALSE;      
   }
 
-  AliTRDcalibDB  *calibration = AliTRDcalibDB::Instance();
+  AliTRDcalibDB* const calibration = AliTRDcalibDB::Instance();
   if (!calibration) {
     AliFatal("No AliTRDcalibDB instance available\n");
     return kFALSE;  
@@ -764,10 +782,12 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
     return kFALSE;
   }
 
-  fMaxThresh            = fReconstructor->GetRecoParam()->GetClusMaxThresh();
-  fSigThresh            = fReconstructor->GetRecoParam()->GetClusSigThresh();
-  fMinMaxCutSigma       = fReconstructor->GetRecoParam()->GetMinMaxCutSigma();
-  fMinLeftRightCutSigma = fReconstructor->GetRecoParam()->GetMinLeftRightCutSigma();
+  const AliTRDrecoParam *const recoParam = fReconstructor->GetRecoParam();
+
+  fMaxThresh            = recoParam->GetClusMaxThresh();
+  fSigThresh            = recoParam->GetClusSigThresh();
+  fMinMaxCutSigma       = recoParam->GetMinMaxCutSigma();
+  fMinLeftRightCutSigma = recoParam->GetMinLeftRightCutSigma();
 
   Int_t istack  = fIndexes->GetStack();
   fLayer  = fIndexes->GetLayer();
@@ -781,6 +801,8 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
     return kFALSE;
   }
 
+  AliDebug(2, Form("Det[%d] @ Sec[%d] Stk[%d] Ly[%d]", fDet, isector, istack, fLayer));
+
   // TRD space point transformation
   fTransform->SetDetector(det);
 
@@ -792,8 +814,31 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
     AddTrackletsToArray();
 
   fColMax    = fDigits->GetNcol();
-  //Int_t nRowMax    = fDigits->GetNrow();
-  fTimeTotal = fDigits->GetNtime();
+  fTimeTotal = fDigitsManager->GetDigitsParam()->GetNTimeBins(det);
+
+  // Check consistency between OCDB and raw data
+  Int_t nTimeOCDB = calibration->GetNumberOfTimeBinsDCS();
+  if(TestBit(kHLT)){
+    if((nTimeOCDB > -1) && (fTimeTotal != nTimeOCDB)){
+      AliWarning(Form("Number of timebins does not match OCDB value (RAW[%d] OCDB[%d]), using raw value"
+                     ,fTimeTotal,nTimeOCDB));
+    }
+  }else{
+    if(nTimeOCDB == -1){
+      AliWarning("Undefined number of timebins in OCDB, using value from raw data.");
+      if(!fTimeTotal>0){
+       AliError("Number of timebins in raw data is negative, skipping chamber!");
+       return kFALSE;
+      }
+    }else if(nTimeOCDB == -2){
+      AliError("Mixed number of timebins in OCDB, no reconstruction of TRD data!"); 
+      return kFALSE;
+    }else if(fTimeTotal != nTimeOCDB){
+      AliError(Form("Number of timebins in raw data does not match OCDB value (RAW[%d] OCDB[%d]), skipping chamber!"
+                   ,fTimeTotal,nTimeOCDB));
+      return kFALSE;
+    }
+  }
 
   // Detector wise calibration object for the gain factors
   const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
@@ -812,15 +857,14 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
   // Calibration object with the pad status
   fCalPadStatusROC       = calibration->GetPadStatusROC(fDet);
 
-  SetBit(kLUT, fReconstructor->GetRecoParam()->UseLUT());
-  SetBit(kGAUS, fReconstructor->GetRecoParam()->UseGAUS());
-  SetBit(kHLT, fReconstructor->IsHLT());
-  
   firstClusterROC = -1;
   fClusterROC     =  0;
 
+  SetBit(kLUT, recoParam->UseLUT());
+  SetBit(kGAUS, recoParam->UseGAUS());
+
   // Apply the gain and the tail cancelation via digital filter
-  if(fReconstructor->GetRecoParam()->UseTailCancelation()) TailCancelation();
+  if(recoParam->UseTailCancelation()) TailCancelation(recoParam);
 
   MaxStruct curr, last;
   Int_t nMaximas = 0, nCorrupted = 0;
@@ -840,7 +884,7 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
   }
   if(last.row>-1) CreateCluster(last);
 
-  if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){
+  if(recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){
     TTreeSRedirector* fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
     (*fDebugStream) << "MakeClusters"
       << "Detector="   << det
@@ -878,14 +922,16 @@ Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Sh
    ,fCalPadStatusROC->GetStatus(Max.col+1, Max.row)
   };
 
-gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
- Signals[0] = (Short_t)((fDigits->GetData(Max.row, Max.col-1, Max.time) - fBaseline) / gain + 0.5f);
   gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
-  Signals[2] = (Short_t)((fDigits->GetData(Max.row, Max.col+1, Max.time) - fBaseline) / gain  + 0.5f);
+  Signals[0] = (Short_t)((fDigits->GetData(Max.row, Max.col-1, Max.time) - fBaseline) / gain + 0.5f);
+  gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+1,Max.row);
+  Signals[2] = (Short_t)((fDigits->GetData(Max.row, Max.col+1, Max.time) - fBaseline) / gain + 0.5f);
 
   if(!(status[0] | status[1] | status[2])) {//all pads are good
     if ((Signals[2] <= Signals[1]) && (Signals[0] <  Signals[1])) {
       if ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh)) {
+       if(Signals[0]<0)Signals[0]=0;
+       if(Signals[2]<0)Signals[2]=0;
         Float_t  noiseSumThresh = fMinLeftRightCutSigma
           * fCalNoiseDetValue
           * fCalNoiseROC->GetValue(Max.col, Max.row);
@@ -895,6 +941,8 @@ gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
       }
     }
   } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
+    if(Signals[0]<0)Signals[0]=0;
+    if(Signals[2]<0)Signals[2]=0;
     if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] >= fSigThresh) { 
       Signals[2]=0;
       SetPadStatus(status[2], padStatus);
@@ -1196,83 +1244,71 @@ Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *con
 }
 
 //_____________________________________________________________________________
-void AliTRDclusterizer::TailCancelation()
+void AliTRDclusterizer::TailCancelation(const AliTRDrecoParam* const recoParam)
 {
   //
-  // Applies the tail cancelation and gain factors: 
-  // Transform fDigits to fDigits
+  // Applies the tail cancelation
   //
 
   Int_t iRow  = 0;
   Int_t iCol  = 0;
   Int_t iTime = 0;
 
-  Double_t *inADC = new Double_t[fTimeTotal];  // ADC data before tail cancellation
-  Double_t *outADC = new Double_t[fTimeTotal];  // ADC data after tail cancellation
+  Float_t *arr = new Float_t[fTimeTotal];  // temp array containing the ADC signals
 
-  fIndexes->ResetCounters();
   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
+  Bool_t debugStreaming = recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming();
+  Int_t nexp = recoParam->GetTCnexp();
   while(fIndexes->NextRCIndex(iRow, iCol))
     {
-      Bool_t corrupted = kFALSE;
-      for (iTime = 0; iTime < fTimeTotal; iTime++) 
-        {        
-          // Apply gain gain factor
-          inADC[iTime]   = fDigits->GetData(iRow,iCol,iTime);
-          if (fCalPadStatusROC->GetStatus(iCol, iRow)) corrupted = kTRUE;
-          outADC[iTime]  = inADC[iTime];
-         if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming()){
-           (*fDebugStream) << "TailCancellation"
-                             << "col="  << iCol
-                             << "row="  << iRow
-                             << "time=" << iTime
-                             << "inADC=" << inADC[iTime]
-                             << "outADC=" << outADC[iTime]
-                             << "corrupted=" << corrupted
-                             << "\n";
-         }
-        }
-      if (!corrupted)
-        {
-          // Apply the tail cancelation via the digital filter
-          // (only for non-coorupted pads)
-         DeConvExp(&inADC[0],&outADC[0],fTimeTotal,fReconstructor->GetRecoParam() ->GetTCnexp());
-        }
-
-      for(iTime = 0; iTime < fTimeTotal; iTime++)//while (fIndexes->NextTbinIndex(iTime))
-        {
-          // Store the amplitude of the digit if above threshold
-          if (outADC[iTime] > 0)
-           fDigits->SetData(iRow,iCol,iTime,TMath::Nint(outADC[iTime]));
-         else
-           fDigits->SetData(iRow,iCol,iTime,0);
-        } // while itime
+      // if corrupted then don't make the tail cancallation
+      if (fCalPadStatusROC->GetStatus(iCol, iRow)) continue;
 
+      // Save data into the temporary processing array and substract the baseline,
+      // since DeConvExp does not expect a baseline
+      for (iTime = 0; iTime < fTimeTotal; iTime++) 
+       arr[iTime] = fDigits->GetData(iRow,iCol,iTime)-fBaseline;
+          
+      if(debugStreaming){
+       for (iTime = 0; iTime < fTimeTotal; iTime++) 
+         (*fDebugStream) << "TailCancellation"
+                         << "col="  << iCol
+                         << "row="  << iRow
+                         << "time=" << iTime
+                         << "arr=" << arr[iTime]
+                         << "\n";
+      }
+      
+      // Apply the tail cancelation via the digital filter
+      DeConvExp(arr,fTimeTotal,nexp);
+
+      // Save tailcancalled data and add the baseline
+      for(iTime = 0; iTime < fTimeTotal; iTime++)
+       fDigits->SetData(iRow,iCol,iTime,(Short_t)(arr[iTime] + fBaseline + 0.5f));
+      
     } // while irow icol
 
-  delete [] inADC;
-  delete [] outADC;
+  delete [] arr;
 
   return;
 
 }
 
 //_____________________________________________________________________________
-void AliTRDclusterizer::DeConvExp(const Double_t *const source, Double_t *const target
-                                 ,const Int_t n, const Int_t nexp) 
+void AliTRDclusterizer::DeConvExp(Float_t *const arr, const Int_t nTime, const Int_t nexp) 
 {
   //
   // Tail cancellation by deconvolution for PASA v4 TRF
   //
 
-  Double_t rates[2];
-  Double_t coefficients[2];
+  Float_t rates[2];
+  Float_t coefficients[2];
 
   // Initialization (coefficient = alpha, rates = lambda)
-  Double_t r1 = 1.0;
-  Double_t r2 = 1.0;
-  Double_t c1 = 0.5;
-  Double_t c2 = 0.5;
+  Float_t r1 = 1.0;
+  Float_t r2 = 1.0;
+  Float_t c1 = 0.5;
+  Float_t c2 = 0.5;
 
   if (nexp == 1) {   // 1 Exponentials
     r1 = 1.156;
@@ -1300,19 +1336,19 @@ void AliTRDclusterizer::DeConvExp(const Double_t *const source, Double_t *const
   Int_t i = 0;
   Int_t k = 0;
 
-  Double_t reminder[2];
-  Double_t correction = 0.0;
-  Double_t result     = 0.0;
+  Float_t reminder[2];
+  Float_t correction = 0.0;
+  Float_t result     = 0.0;
 
   // Attention: computation order is important
   for (k = 0; k < nexp; k++) {
     reminder[k] = 0.0;
   }
 
-  for (i = 0; i < n; i++) {
+  for (i = 0; i < nTime; i++) {
 
-    result    = (source[i] - correction);    // No rescaling
-    target[i] = result;
+    result = (arr[i] - correction);    // No rescaling
+    arr[i] = result;
 
     for (k = 0; k < nexp; k++) {
       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
@@ -1335,8 +1371,9 @@ void AliTRDclusterizer::ResetRecPoints()
   //
 
   if (fRecPoints) {
-    fRecPoints->Delete();
-    delete fRecPoints;
+    fRecPoints->Clear();
+    fNoOfClusters = 0;
+    //    delete fRecPoints;
   }
 }