]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDclusterizer.cxx
Moveing AliNeutralTrackParam to STEERBase to avoid dependence of libAOD on libESD
[u/mrichter/AliRoot.git] / TRD / AliTRDclusterizer.cxx
index 1d215d4968b7034088e8b931cd5377d6a86f06e4..6611a8276d9c829a9e19b4914aa1d6276f506598 100644 (file)
@@ -26,6 +26,7 @@
 
 #include "AliRunLoader.h"
 #include "AliLoader.h"
+#include "AliTreeLoader.h"
 #include "AliAlignObj.h"
 
 #include "AliTRDclusterizer.h"
@@ -68,6 +69,7 @@ AliTRDclusterizer::AliTRDclusterizer(const AliTRDReconstructor *const rec)
   ,fDigits(NULL)
   ,fIndexes(NULL)
   ,fMaxThresh(0)
+  ,fMaxThreshTest(0)
   ,fSigThresh(0)
   ,fMinMaxCutSigma(0)
   ,fMinLeftRightCutSigma(0)
@@ -128,6 +130,7 @@ AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title, co
   ,fDigits(NULL)
   ,fIndexes(NULL)
   ,fMaxThresh(0)
+  ,fMaxThreshTest(0)
   ,fSigThresh(0)
   ,fMinMaxCutSigma(0)
   ,fMinLeftRightCutSigma(0)
@@ -183,6 +186,7 @@ AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c)
   ,fDigits(NULL)
   ,fIndexes(NULL)
   ,fMaxThresh(0)
+  ,fMaxThreshTest(0)
   ,fSigThresh(0)
   ,fMinMaxCutSigma(0)
   ,fMinLeftRightCutSigma(0)
@@ -287,6 +291,7 @@ void AliTRDclusterizer::Copy(TObject &c) const
   ((AliTRDclusterizer &) c).fDigits      = NULL;
   ((AliTRDclusterizer &) c).fIndexes       = NULL;
   ((AliTRDclusterizer &) c).fMaxThresh     = 0;
+  ((AliTRDclusterizer &) c).fMaxThreshTest = 0;
   ((AliTRDclusterizer &) c).fSigThresh     = 0;
   ((AliTRDclusterizer &) c).fMinMaxCutSigma= 0;
   ((AliTRDclusterizer &) c).fMinLeftRightCutSigma = 0;
@@ -371,46 +376,6 @@ Bool_t AliTRDclusterizer::OpenOutput(TTree *const clusterTree)
   return kTRUE;
 }
 
-//_____________________________________________________________________________
-Bool_t AliTRDclusterizer::OpenTrackletOutput()
-{
-  //
-  // Tracklet writing
-  //
-
-  if (fReconstructor->IsWritingTracklets()){
-    TString evfoldname = AliConfig::GetDefaultEventFolderName();
-    fRunLoader         = AliRunLoader::GetRunLoader(evfoldname);
-
-    if (!fRunLoader) {
-      fRunLoader = AliRunLoader::Open("galice.root");
-    }
-    if (!fRunLoader) {
-      AliError(Form("Can not open session for file galice.root."));
-      return kFALSE;
-    }
-
-    UInt_t **leaves = new UInt_t *[2];
-    AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
-    if (!dl) {
-      AliError("Could not get the tracklets data loader!");
-      dl = new AliDataLoader("TRD.Tracklets.root","tracklets", "tracklets");
-      fRunLoader->GetLoader("TRDLoader")->AddDataLoader(dl);
-    }
-    fTrackletTree = dl->Tree();
-    if (!fTrackletTree)
-      {
-       dl->MakeTree();
-       fTrackletTree = dl->Tree();
-      }
-    TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
-    if (!trkbranch)
-      fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
-  }
-
-  return kTRUE;
-}
-
 //_____________________________________________________________________________
 Bool_t AliTRDclusterizer::OpenInput(Int_t nEvent)
 {
@@ -511,8 +476,8 @@ Bool_t AliTRDclusterizer::WriteTracklets(Int_t det)
     }
   }
 
-  AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
-  dl->WriteData("OVERWRITE");
+//  AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets"); //jkl: wrong here
+//  dl->WriteData("OVERWRITE"); //jkl: wrong here
   //dl->Unload();
   delete [] leaves;
 
@@ -611,12 +576,9 @@ 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);
   }
-  fReconstructor->SetDigitsParam(fDigitsManager->GetDigitsParam());
   
   if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
 
@@ -653,6 +615,23 @@ Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
 
   fDigitsManager->SetUseDictionaries(TestBit(kLabels));
 
+  // ----- preparing tracklet output -----
+  if (fReconstructor->IsWritingTracklets()) {
+    AliDataLoader *trklLoader = AliRunLoader::Instance()->GetLoader("TRDLoader")->GetDataLoader("tracklets");
+    if (!trklLoader) {
+      //AliInfo("Could not get the tracklets data loader, adding it now!");
+      trklLoader = new AliDataLoader("TRD.Tracklets.root","tracklets", "tracklets");
+      AliRunLoader::Instance()->GetLoader("TRDLoader")->AddDataLoader(trklLoader);
+    }
+    AliTreeLoader *trklTreeLoader = dynamic_cast<AliTreeLoader*> (trklLoader->GetBaseLoader("tracklets-raw"));
+    if (!trklTreeLoader) {
+      trklTreeLoader = new AliTreeLoader("tracklets-raw", trklLoader);
+      trklLoader->AddBaseLoader(trklTreeLoader);
+    }
+    if (!trklTreeLoader->Tree())
+      trklTreeLoader->MakeTree();
+  }
+
   // tracklet container for raw tracklet writing
   if (!fTrackletContainer && ( fReconstructor->IsWritingTracklets() || fReconstructor->IsProcessingTracklets() )) {
     // maximum tracklets for one HC
@@ -660,6 +639,8 @@ Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
     fTrackletContainer = new UInt_t *[2];
     fTrackletContainer[0] = new UInt_t[kTrackletChmb]; 
     fTrackletContainer[1] = new UInt_t[kTrackletChmb]; 
+    memset(fTrackletContainer[0], 0, kTrackletChmb*sizeof(UInt_t)); //jkl
+    memset(fTrackletContainer[1], 0, kTrackletChmb*sizeof(UInt_t)); //jkl
   }
 
   if(!fRawStream)
@@ -686,8 +667,17 @@ Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
     if (!fReconstructor->IsWritingTracklets()) continue;
     if (*(fTrackletContainer[0]) > 0 || *(fTrackletContainer[1]) > 0) WriteTracklets(det);
   }
-  fReconstructor->SetDigitsParam(fDigitsManager->GetDigitsParam());
   
+  if (fReconstructor->IsWritingTracklets()) {
+    if (AliDataLoader *trklLoader = AliRunLoader::Instance()->GetLoader("TRDLoader")->GetDataLoader("tracklets")) {
+      if (trklLoader) {
+       if (AliTreeLoader *trklTreeLoader = (AliTreeLoader*) trklLoader->GetBaseLoader("tracklets-raw"))
+         trklTreeLoader->WriteData("OVERWRITE");
+       trklLoader->UnloadAll();
+      }
+    }
+  }
+
   if (fTrackletContainer){
     delete [] fTrackletContainer[0];
     delete [] fTrackletContainer[1];
@@ -764,7 +754,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;
@@ -788,10 +778,12 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
 
   const AliTRDrecoParam *const recoParam = fReconstructor->GetRecoParam();
 
-  fMaxThresh            = recoParam->GetClusMaxThresh();
-  fSigThresh            = recoParam->GetClusSigThresh();
+  fMaxThresh            = (Short_t)recoParam->GetClusMaxThresh();
+  fMaxThreshTest        = (Short_t)(recoParam->GetClusMaxThresh()/2+fBaseline);
+  fSigThresh            = (Short_t)recoParam->GetClusSigThresh();
   fMinMaxCutSigma       = recoParam->GetMinMaxCutSigma();
   fMinLeftRightCutSigma = recoParam->GetMinLeftRightCutSigma();
+  const Int_t iEveryNTB = recoParam->GetRecEveryNTB();
 
   Int_t istack  = fIndexes->GetStack();
   fLayer  = fIndexes->GetLayer();
@@ -818,8 +810,7 @@ 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();
@@ -876,9 +867,9 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
 
   // Here the clusterfining is happening
   
-  for(curr.time = 0; curr.time < fTimeTotal; curr.time++){
+  for(curr.time = 0; curr.time < fTimeTotal; curr.time+=iEveryNTB){
     while(fIndexes->NextRCIndex(curr.row, curr.col)){
-      if(IsMaximum(curr, curr.padStatus, &curr.signals[0])){
+      if(fDigits->GetData(curr.row, curr.col, curr.time) > fMaxThreshTest && IsMaximum(curr, curr.padStatus, &curr.signals[0])){
         if(last.row>-1){
           if(curr.time==last.time && curr.row==last.row && curr.col==last.col+2) FivePadCluster(last, curr);
           CreateCluster(last);
@@ -914,10 +905,10 @@ Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Sh
 
   Float_t gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col,Max.row);
   Signals[1] = (Short_t)((fDigits->GetData(Max.row, Max.col, Max.time) - fBaseline) / gain + 0.5f);
-  if(Signals[1] < fMaxThresh) return kFALSE;
+  if(Signals[1] <= fMaxThresh) return kFALSE;
 
-  Float_t  noiseMiddleThresh = fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.col, Max.row);
-  if (Signals[1] < noiseMiddleThresh) return kFALSE;
+  Short_t noiseMiddleThresh = (Short_t)(fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.col, Max.row));
+  if (Signals[1] <= noiseMiddleThresh) return kFALSE;
 
   if (Max.col + 1 >= fColMax || Max.col < 1) return kFALSE;
 
@@ -934,13 +925,12 @@ Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Sh
 
   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[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);
-        if ((Signals[2]+Signals[0]+Signals[1]) < noiseSumThresh) return kFALSE;
+        Short_t noiseSumThresh = (Short_t)(fMinLeftRightCutSigma * fCalNoiseDetValue
+                                          * fCalNoiseROC->GetValue(Max.col, Max.row));
+        if ((Signals[2]+Signals[0]+Signals[1]) <= noiseSumThresh) return kFALSE;
         padStatus = 0;
         return kTRUE;
       }
@@ -948,18 +938,18 @@ Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Sh
   } 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) { 
+    if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] > fSigThresh) { 
       Signals[2]=0;
       SetPadStatus(status[2], padStatus);
       return kTRUE;
     } 
-    else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] >= fSigThresh) {
+    else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] > fSigThresh) {
       Signals[0]=0;
       SetPadStatus(status[0], padStatus);
       return kTRUE;
     }
-    else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh))) {
-      Signals[1] = (Short_t)(fMaxThresh + 0.5f);
+    else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh))) {
+      Signals[1] = fMaxThresh;
       SetPadStatus(status[1], padStatus);
       return kTRUE;
     }
@@ -1039,7 +1029,7 @@ void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
   cluster.SetLabel(Max.time,2);
 
   //needed for HLT reconstruction
-  AddClusterToArray(&cluster); 
+  AddClusterToArray(&cluster);
 
   // Store the index of the first cluster in the current ROC
   if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
@@ -1259,8 +1249,6 @@ void AliTRDclusterizer::TailCancelation(const AliTRDrecoParam* const recoParam)
   Int_t iCol  = 0;
   Int_t iTime = 0;
 
-  Float_t *arr = new Float_t[fTimeTotal];  // temp array containing the ADC signals
-
   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
   Bool_t debugStreaming = recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming();
   Int_t nexp = recoParam->GetTCnexp();
@@ -1269,38 +1257,25 @@ void AliTRDclusterizer::TailCancelation(const AliTRDrecoParam* const recoParam)
       // 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";
+       for (iTime = 0; iTime < fTimeTotal; iTime++) 
+         (*fDebugStream) << "TailCancellation"
+                         << "col="  << iCol
+                         << "row="  << iRow
+                         << "\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));
+      DeConvExp(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp);
       
     } // while irow icol
 
-  delete [] arr;
-
   return;
 
 }
 
 //_____________________________________________________________________________
-void AliTRDclusterizer::DeConvExp(Float_t *const arr, const Int_t nTime, const Int_t nexp) 
+void AliTRDclusterizer::DeConvExp(Short_t *const arr, const Int_t nTime, const Int_t nexp)
 {
   //
   // Tail cancellation by deconvolution for PASA v4 TRF
@@ -1352,8 +1327,8 @@ void AliTRDclusterizer::DeConvExp(Float_t *const arr, const Int_t nTime, const I
 
   for (i = 0; i < nTime; i++) {
 
-    result = (arr[i] - correction);    // No rescaling
-    arr[i] = result;
+    result = arr[i] - correction - fBaseline;    // No rescaling
+    arr[i] = (Short_t)(result + fBaseline + 0.5f);
 
     for (k = 0; k < nexp; k++) {
       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
@@ -1376,8 +1351,9 @@ void AliTRDclusterizer::ResetRecPoints()
   //
 
   if (fRecPoints) {
-    fRecPoints->Delete();
-    delete fRecPoints;
+    fRecPoints->Clear();
+    fNoOfClusters = 0;
+    //    delete fRecPoints;
   }
 }