changes in AliTRDdigitsManager
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 23 Jan 2009 11:58:35 +0000 (11:58 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 23 Jan 2009 11:58:35 +0000 (11:58 +0000)
- important correction in the destructor
- TObjArrays may contain only one entry if raw data is reconstructed. faster!
- cleaned up code

changes in AliTRDclusterizer
- important correction concering signals[7] in CreateCluster
- some corrections of the padStatus handeling in IsMaximum
- make use of changes in digitsManager
- cleaned up code of IsMaximum
- preparation of the code for upcoming SIMD implementation

changes in AliTRDReconstructor
- make use of changes in digitsManager

TRD/AliTRDReconstructor.cxx
TRD/AliTRDclusterizer.cxx
TRD/AliTRDclusterizer.h
TRD/AliTRDdigitsManager.cxx
TRD/AliTRDdigitsManager.h

index cc735de..b1e8da1 100644 (file)
@@ -157,9 +157,9 @@ void AliTRDReconstructor::Reconstruct(AliRawReader *rawReader
   rawReader->Select("TRD");
 
   // New (fast) cluster finder
-  AliTRDclusterizer clusterer("clusterer","TRD clusterizer");
-  clusterer.SetReconstructor(this);
-  clusterer.OpenOutput(clusterTree);
+  AliTRDclusterizer clusterer("clusterer","TRD clusterizer",this);
+  //clusterer.SetReconstructor(this);                 //     ^| "this" tells the digitsmanager that we are reading raw files
+  clusterer.OpenOutput(clusterTree);                  //        it is not strictly necessaray but will give a speed up
   clusterer.SetAddLabels(kFALSE);
   clusterer.Raw2ClustersChamber(rawReader);
   
@@ -181,8 +181,8 @@ void AliTRDReconstructor::Reconstruct(TTree *digitsTree
   //AliInfo("Reconstruct TRD clusters from Digits [Digit TTree -> Cluster TTree]");
 
   AliTRDclusterizer clusterer("clusterer","TRD clusterizer");
-  clusterer.SetReconstructor(this);
-  clusterer.OpenOutput(clusterTree);
+  clusterer.SetReconstructor(this);                  //    ^| no this, because we are reading from digitsTree
+  clusterer.OpenOutput(clusterTree);                 //       it is necessary to NOT have the "this" here!
   clusterer.ReadDigits(digitsTree);
   clusterer.MakeClusters();
 
index 85f2cbd..d904886 100644 (file)
 ClassImp(AliTRDclusterizer)
 
 //_____________________________________________________________________________
-AliTRDclusterizer::AliTRDclusterizer(AliTRDReconstructor *rec)
+AliTRDclusterizer::AliTRDclusterizer(const AliTRDReconstructor *const rec)
   :TNamed()
   ,fReconstructor(rec)  
   ,fRunLoader(NULL)
   ,fClusterTree(NULL)
   ,fRecPoints(NULL)
   ,fTrackletTree(NULL)
-  ,fDigitsManager(NULL)
+  ,fDigitsManager(new AliTRDdigitsManager(rec))
   ,fTrackletContainer(NULL)
   ,fAddLabels(kTRUE)
   ,fRawVersion(2)
@@ -117,14 +117,14 @@ AliTRDclusterizer::AliTRDclusterizer(AliTRDReconstructor *rec)
 }
 
 //_____________________________________________________________________________
-AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title, AliTRDReconstructor *rec)
+AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title, const AliTRDReconstructor *const rec)
   :TNamed(name,title)
   ,fReconstructor(rec)
   ,fRunLoader(NULL)
   ,fClusterTree(NULL)
   ,fRecPoints(NULL)
   ,fTrackletTree(NULL)
-  ,fDigitsManager(new AliTRDdigitsManager())
+  ,fDigitsManager(new AliTRDdigitsManager(rec))
   ,fTrackletContainer(NULL)
   ,fAddLabels(kTRUE)
   ,fRawVersion(2)
@@ -631,7 +631,7 @@ Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
 
   // Create the digits manager
   if (!fDigitsManager){
-    fDigitsManager = new AliTRDdigitsManager();
+    fDigitsManager = new AliTRDdigitsManager(fReconstructor);
     fDigitsManager->CreateArrays();
   }
 
@@ -646,7 +646,6 @@ Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
     fTrackletContainer[1] = new UInt_t[kTrackletChmb]; 
   }
 
-
   AliTRDrawStreamBase *input = AliTRDrawStreamBase::GetRawStream(rawReader);
 
   AliInfo(Form("Stream version: %s", input->IsA()->GetName()));
@@ -655,16 +654,15 @@ Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
   while ((det = input->NextChamber(fDigitsManager,fTrackletContainer)) >= 0){
     Bool_t iclusterBranch = kFALSE;
     if (fDigitsManager->GetIndexes(det)->HasEntry()){
-    iclusterBranch = MakeClusters(det);
+      iclusterBranch = MakeClusters(det);
     }
 
-    fDigitsManager->RemoveDigits(det);
-    fDigitsManager->RemoveDictionaries(det);      
-    fDigitsManager->ClearIndexes(det);
-
+    fDigitsManager->ResetArrays(det);
+    
     if (!fReconstructor->IsWritingTracklets()) continue;
     if (*(fTrackletContainer[0]) > 0 || *(fTrackletContainer[1]) > 0) WriteTracklets(det);
   }
+  
   if (fReconstructor->IsWritingTracklets()){
     delete [] fTrackletContainer[0];
     delete [] fTrackletContainer[1];
@@ -823,36 +821,36 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
   // Apply the gain and the tail cancelation via digital filter
   TailCancelation();   
 
-  ClusterizerStruct curr, last;
+  MaxStruct curr, last;
   last.Row = -1;
   Int_t nMaximas = 0, nCorrupted = 0;
-  Double_t Ratio = 1;
+  Float_t Ratio = 1;
 
   // Here the clusterfining is happening
   
   for(curr.Time = 0; curr.Time < fTimeTotal; curr.Time++)
     while(fIndexes->NextRCIndex(curr.Row, curr.Col))
-      if(IsMaximum(curr.Row, curr.Col, curr.Time, curr.padStatus, &curr.Signals[0]))
+      if(IsMaximum(curr, curr.padStatus, &curr.Signals[0]))
        {
          if(last.Row>-1)
+           {
+             last.Signals[0] *= Ratio;
+             if(curr.Row==last.Row && curr.Col==last.Col+2)
                {
-                 last.Signals[0] *= Ratio;
-                 if(curr.Row==last.Row && curr.Col==last.Col+2)
+                 if(IsFivePadCluster(last, curr, Ratio))
                    {
-                     if(IsFivePadCluster(last.Row, last.Col, last.Time, &last.Signals[0], &curr.Signals[0], Ratio))
-                       {
-                         last.Signals[2] *= Ratio;
-                         Ratio = 1 - Ratio;
-                       }else Ratio = 1;
+                     last.Signals[2] *= Ratio;
+                     Ratio = 1 - Ratio;
                    }else Ratio = 1;
-                 CreateCluster(last.Row, last.Col, last.Time, &last.Signals[0], last.padStatus);
-               }
+               }else Ratio = 1;
+             CreateCluster(last);
+           }
          last=curr;
        }
   if(last.Row>-1)
     {
       last.Signals[0] *= Ratio;
-      CreateCluster(last.Row, last.Col, last.Time, &last.Signals[0], last.padStatus);
+      CreateCluster(last);
     }
 
   if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kClusterizer) > 2){
@@ -873,66 +871,54 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
 }
 
 //_____________________________________________________________________________
-Bool_t AliTRDclusterizer::IsMaximum(const Int_t row, const Int_t col, const Int_t time, 
-                                   UChar_t &pasStatus, Double_t *const Signals) 
+Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Float_t *const Signals) 
 {
   //
   // Returns true if this row,col,time combination is a maximum. 
   // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
   //
 
-  Signals[1] = fDigitsOut->GetData(row,col,time);
+  Signals[1] = fDigitsOut->GetData(Max.Row, Max.Col, Max.Time);
   if(Signals[1] < fMaxThresh) return kFALSE;
 
-  Float_t  noiseMiddleThresh = fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(col,row);
+  Float_t  noiseMiddleThresh = fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.Col, Max.Row);
   if (Signals[1] < noiseMiddleThresh) return kFALSE;
 
-  if (col + 1 >= fColMax || col < 1) return kFALSE;
-  UChar_t status[3]={0, 0, 0};
-  pasStatus = 0;
+  if (Max.Col + 1 >= fColMax || Max.Col < 1) return kFALSE;
 
-  status[1] = fDigitsIn->GetPadStatus(row,col,time);
-  //if(status[1]) SETBIT(pasStatus, AliTRDcluster::kMaskedCenter);//TR: mod: this is already done by SetPadStatus
+  UChar_t status[3]={fDigitsIn->GetPadStatus(Max.Row, Max.Col-1, Max.Time), 
+                    fDigitsIn->GetPadStatus(Max.Row, Max.Col,   Max.Time), 
+                    fDigitsIn->GetPadStatus(Max.Row, Max.Col+1, Max.Time)};
 
-  Signals[2] = fDigitsOut->GetData(row,col+1,time);
-  status[2] = fDigitsIn->GetPadStatus(row,col+1,time);
-  //if(status[2]) SETBIT(pasStatus, AliTRDcluster::kMaskedLeft);//TR: mod: this is already done by SetPadStatus
-    
-  Signals[0] = fDigitsOut->GetData(row,col-1,time);
-  status[0] = fDigitsIn->GetPadStatus(row,col-1,time);
-  //if(status[0]) SETBIT(pasStatus, AliTRDcluster::kMaskedRight);//TR: mod: this is already done by SetPadStatus
-    
-  // reject candidates with more than 1 problematic pad
-  if(pasStatus >= 3) return kFALSE;
-    
-  if (!status[1]) { // good central pad
-    if (!pasStatus) { // all pads are OK
-      if ((Signals[2] <= Signals[1]) && (Signals[0] <  Signals[1])) {
-       if ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh)) {
-         Float_t  noiseSumThresh = fMinLeftRightCutSigma
-           * fCalNoiseDetValue
-           * fCalNoiseROC->GetValue(col,row);
-         if ((Signals[2]+Signals[0]+Signals[1]) >= noiseSumThresh)
-           return kTRUE;
-       }
-      }
-    } else { // one of the neighbouring pads are bad
-      if (status[2] && Signals[0] < Signals[1] && Signals[0] >= fSigThresh) { 
-       fDigitsOut->SetData(row, col+1, time, 0.);//TR: mod: was: SetData(row, col, time+1, 0.)
-       SetPadStatus(status[2], pasStatus);
-       return kTRUE;
-      } 
-      else if (status[0] && Signals[2] <= Signals[1] && Signals[2] >= fSigThresh) { 
-       fDigitsOut->SetData(row, col-1, time, 0.);//TR: mod: was: SetData(row, col, time-1, 0.)
-       SetPadStatus(status[0], pasStatus);
+  Signals[0] = fDigitsOut->GetData(Max.Row, Max.Col-1, Max.Time);
+  Signals[2] = fDigitsOut->GetData(Max.Row, Max.Col+1, Max.Time);  
+
+  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)) {
+       Float_t  noiseSumThresh = fMinLeftRightCutSigma
+         * fCalNoiseDetValue
+         * fCalNoiseROC->GetValue(Max.Col, Max.Row);
+       if ((Signals[2]+Signals[0]+Signals[1]) < noiseSumThresh) return kFALSE;
+       padStatus = 0;
        return kTRUE;
       }
     }
-  } 
-  else { // wrong maximum pad
-    if ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh)) {
-      fDigitsOut->SetData(row,col,time,fMaxThresh);
-      SetPadStatus(status[1], pasStatus);
+  }
+  else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
+    if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] >= fSigThresh) { 
+      fDigitsOut->SetData(Max.Row, Max.Col+1, Max.Time, 0.);
+      SetPadStatus(status[2], padStatus);
+      return kTRUE;
+    } 
+    else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] >= fSigThresh) {
+      fDigitsOut->SetData(Max.Row, Max.Col-1, Max.Time, 0.);
+      SetPadStatus(status[0], padStatus);
+      return kTRUE;
+    }
+    else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh))) {
+      fDigitsOut->SetData(Max.Row, Max.Col, Max.Time, fMaxThresh);
+      SetPadStatus(status[1], padStatus);
       return kTRUE;
     }
   }
@@ -940,41 +926,37 @@ Bool_t AliTRDclusterizer::IsMaximum(const Int_t row, const Int_t col, const Int_
 }
 
 //_____________________________________________________________________________
-Bool_t AliTRDclusterizer::IsFivePadCluster(const Int_t row, const Int_t col, const Int_t time, 
-                                          Double_t *SignalsThisMax, Double_t *SignalsNeighbourMax, Double_t &ratio)
+Bool_t AliTRDclusterizer::IsFivePadCluster(const MaxStruct &ThisMax, const MaxStruct &NeighbourMax, Float_t &ratio)
 {
   //
   // Look for 5 pad cluster with minimum in the middle
   // Gives back the ratio
   //
 
-  if (col < fColMax - 3){
-    if (col < fColMax - 5){
-      if (fDigitsOut->GetData(row,col+4,time) >= fSigThresh)
-       return kFALSE;
-    }
-    if (col > 1) {
-      if (fDigitsOut->GetData(row,col-2,time) >= fSigThresh)
-       return kFALSE;
-       }
-    
-  //if (fSignalsThisMax[1] >= 0){ //TR: mod
-
-    const Float_t kEpsilon = 0.01;
-    Double_t padSignal[5] = {SignalsThisMax[0],SignalsThisMax[1],SignalsThisMax[2],
-                            SignalsNeighbourMax[1], SignalsNeighbourMax[2]};
-    
-    // Unfold the two maxima and set the signal on 
-    // the overlapping pad to the ratio
-    ratio = Unfold(kEpsilon,fLayer,padSignal);
-    return kTRUE;
+  if (ThisMax.Col >= fColMax - 3) return kFALSE;
+  if (ThisMax.Col < fColMax - 5){
+    if (fDigitsOut->GetData(ThisMax.Row, ThisMax.Col+4, ThisMax.Time) >= fSigThresh)
+      return kFALSE;
   }
-  return kFALSE;
+  if (ThisMax.Col > 1) {
+    if (fDigitsOut->GetData(ThisMax.Row, ThisMax.Col-2, ThisMax.Time) >= fSigThresh)
+      return kFALSE;
+  }
+  
+  //if (fSignalsThisMax[1] >= 0){ //TR: mod
+  
+  const Float_t kEpsilon = 0.01;
+  Double_t padSignal[5] = {ThisMax.Signals[0], ThisMax.Signals[1], ThisMax.Signals[2],
+                          NeighbourMax.Signals[1], NeighbourMax.Signals[2]};
+  
+  // Unfold the two maxima and set the signal on 
+  // the overlapping pad to the ratio
+  ratio = Unfold(kEpsilon,fLayer,padSignal);
+  return kTRUE;
 }
 
 //_____________________________________________________________________________
-void AliTRDclusterizer::CreateCluster(const Int_t row, const Int_t col, const Int_t time, 
-                                       const Double_t* const clusterSignal, const UChar_t pasStatus)
+void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
 {
   //
   // Creates a cluster at the given position and saves it in fRecPoint
@@ -988,23 +970,23 @@ void AliTRDclusterizer::CreateCluster(const Int_t row, const Int_t col, const In
   if (fReconstructor->GetRecoParam()->IsLUT()) {
     // Calculate the position of the cluster by using the
     // lookup table method
-    clusterPosCol = LUTposition(fLayer,clusterSignal[0]
-                               ,clusterSignal[1]
-                               ,clusterSignal[2]);
+    clusterPosCol = LUTposition(fLayer,Max.Signals[0]
+                               ,Max.Signals[1]
+                               ,Max.Signals[2]);
   } 
   else {
     // Calculate the position of the cluster by using the
     // center of gravity method
-    padSignal[1] = clusterSignal[0];
-    padSignal[2] = clusterSignal[1];
-    padSignal[3] = clusterSignal[2];
-    if(col > 2){
-      padSignal[0] = fDigitsOut->GetData(row,col-2,time);
+    padSignal[1] = Max.Signals[0];
+    padSignal[2] = Max.Signals[1];
+    padSignal[3] = Max.Signals[2];
+    if(Max.Col > 2){
+      padSignal[0] = fDigitsOut->GetData(Max.Row, Max.Col-2, Max.Time);
       if(padSignal[0]>= padSignal[1])
        padSignal[0] = 0;
     }
-    if(col < fColMax - 3){
-      padSignal[4] = fDigitsOut->GetData(row,col+2,time);
+    if(Max.Col < fColMax - 3){
+      padSignal[4] = fDigitsOut->GetData(Max.Row, Max.Col+2, Max.Time);
       if(padSignal[4]>= padSignal[3])
        padSignal[4] = 0;
     }
@@ -1015,34 +997,34 @@ void AliTRDclusterizer::CreateCluster(const Int_t row, const Int_t col, const In
   Int_t nPadCount = 1;
   // Look to the right
   Int_t ii = 1;
-  while (fDigitsOut->GetData(row, col-ii, time) >= fSigThresh) {
+  while (fDigitsOut->GetData(Max.Row, Max.Col-ii, Max.Time) >= fSigThresh) {
     nPadCount++;
     ii++;
-    if (col-ii < 0) break;
+    if (Max.Col < ii) break;
   }
   // Look to the left
   ii = 1;
-  while (fDigitsOut->GetData(row, col+ii, time) >= fSigThresh) {
+  while (fDigitsOut->GetData(Max.Row, Max.Col+ii, Max.Time) >= fSigThresh) {
     nPadCount++;
     ii++;
-    if (col+ii >= fColMax) break;
+    if (Max.Col+ii >= fColMax) break;
   }
 
   // Store the amplitudes of the pads in the cluster for later analysis
   // and check whether one of these pads is masked in the database
-  Short_t signals[7] = { 0, 0, 0, 0, 0, 0, 0 };
-  for(Int_t i = 0; i<3; i++)
-    signals[i+2] = TMath::Nint(clusterSignal[i]);
+  Short_t signals[7] = { 0, 0, TMath::Nint(Max.Signals[0]),
+                              TMath::Nint(Max.Signals[1]), 
+                              TMath::Nint(Max.Signals[2]), 0, 0 };
   for(Int_t i = 0; i<2; i++)
     {
-      if(col+i >= 3)
-       signals[i] = TMath::Nint(fDigitsOut->GetData(row,col-3+i,time));
-      if(col+3-i < fColMax)
-       signals[7-i] = TMath::Nint(fDigitsOut->GetData(row,col+3-i,time));
+      if(Max.Col+i >= 3)
+       signals[i] = TMath::Nint(fDigitsOut->GetData(Max.Row, Max.Col-3+i, Max.Time));
+      if(Max.Col+3-i < fColMax)
+       signals[6-i] = TMath::Nint(fDigitsOut->GetData(Max.Row, Max.Col+3-i, Max.Time));
     }
-  /*for (Int_t jPad = col-3; jPad <= col+3; jPad++) {
+  /*for (Int_t jPad = Max.Col-3; jPad <= Max.Col+3; jPad++) {
     if ((jPad >= 0) && (jPad < fColMax))
-      signals[jPad-col+3] = TMath::Nint(fDigitsOut->GetData(row,jPad,time));
+      signals[jPad-Max.Col+3] = TMath::Nint(fDigitsOut->GetData(Max.Row,jPad,Max.Time));
       }*/
 
   // Transform the local cluster coordinates into calibrated 
@@ -1050,18 +1032,18 @@ void AliTRDclusterizer::CreateCluster(const Int_t row, const Int_t col, const In
   // Here the calibration for T0, Vdrift and ExB is applied as well.
   Double_t clusterXYZ[6];
   clusterXYZ[0] = clusterPosCol;
-  clusterXYZ[1] = clusterSignal[2];
-  clusterXYZ[2] = clusterSignal[1];
-  clusterXYZ[3] = clusterSignal[0];
+  clusterXYZ[1] = Max.Signals[2];
+  clusterXYZ[2] = Max.Signals[1];
+  clusterXYZ[3] = Max.Signals[0];
   clusterXYZ[4] = 0.0;
   clusterXYZ[5] = 0.0;
   Int_t    clusterRCT[3];
-  clusterRCT[0] = row;
-  clusterRCT[1] = col;
+  clusterRCT[0] = Max.Row;
+  clusterRCT[1] = Max.Col;
   clusterRCT[2] = 0;
 
   Bool_t out = kTRUE;
-  if (fTransform->Transform(clusterXYZ,clusterRCT,((UInt_t) time),out,0)) {
+  if (fTransform->Transform(clusterXYZ,clusterRCT,((UInt_t) Max.Time),out,0)) {
 
     // Add the cluster to the output array
     // The track indices will be stored later 
@@ -1082,23 +1064,22 @@ void AliTRDclusterizer::CreateCluster(const Int_t row, const Int_t col, const In
                                                                  0x0,
                                                                  ((Char_t) nPadCount),
                                                                  signals,
-                                                                 ((UChar_t) col), ((UChar_t) row), ((UChar_t) time),
+                                                                 ((UChar_t) Max.Col), ((UChar_t) Max.Row), ((UChar_t) Max.Time),
                                                                  clusterTimeBin, clusterPosCol,
                                                                  fVolid);
     cluster->SetInChamber(!out);
 
-    UChar_t maskPosition = GetCorruption(pasStatus);
-    UChar_t padstatus = GetPadStatus(pasStatus);
+    UChar_t maskPosition = GetCorruption(Max.padStatus);
     if (maskPosition) { 
       cluster->SetPadMaskedPosition(maskPosition);
-      cluster->SetPadMaskedStatus(padstatus);
+      cluster->SetPadMaskedStatus(GetPadStatus(Max.padStatus));
     }
 
-    // Temporarily store the row, column and time bin of the center pad
+    // Temporarily store the Max.Row, column and time bin of the center pad
     // Used to later on assign the track indices
-    cluster->SetLabel(row, 0);
-    cluster->SetLabel(col, 1);
-    cluster->SetLabel(time,2);
+    cluster->SetLabel(Max.Row, 0);
+    cluster->SetLabel(Max.Col, 1);
+    cluster->SetLabel(Max.Time,2);
 
     // Store the index of the first cluster in the current ROC
     if (firstClusterROC < 0) {
@@ -1204,7 +1185,7 @@ Double_t AliTRDclusterizer::GetCOG(Double_t signal[5]) const
 }
 
 //_____________________________________________________________________________
-Double_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, Double_t *padSignal) const
+Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, Double_t *padSignal) const
 {
   //
   // Method to unfold neighbouring maxima.
index b4b7f2b..bb2ca65 100644 (file)
@@ -40,9 +40,23 @@ class AliTRDclusterizer : public TNamed
   enum{
     kOwner = BIT(14)
   };
+
+  struct MaxStruct
+  {
+    Int_t       Row;
+    Int_t       Col;
+    Int_t       Time;
+    UChar_t     padStatus;
+    Float_t     Signals[3];
+    MaxStruct():Row(0),Col(0),Time(0),padStatus(0)
+      {}
+    MaxStruct &operator=(const MaxStruct &a)
+      {Row=a.Row; Col=a.Col; Time=a.Time; padStatus=a.padStatus;
+       memcpy(Signals, a.Signals, 3*sizeof(Signals[0])); return *this;}
+  };
   
-  AliTRDclusterizer(AliTRDReconstructor *rec = 0x0);
-  AliTRDclusterizer(const Text_t* name, const Text_t* title, AliTRDReconstructor *rec = 0x0);
+  AliTRDclusterizer(const AliTRDReconstructor *const rec = 0x0);
+  AliTRDclusterizer(const Text_t* name, const Text_t* title, const AliTRDReconstructor *const rec = 0x0);
   AliTRDclusterizer(const AliTRDclusterizer &c);
   virtual         ~AliTRDclusterizer();
   AliTRDclusterizer &operator=(const AliTRDclusterizer &c);
@@ -84,7 +98,7 @@ class AliTRDclusterizer : public TNamed
                             ,const Int_t nTimeTotal, const Int_t nexp);
   void             TailCancelation();
 
-  virtual Double_t Unfold(Double_t eps, Int_t layer, Double_t *padSignal) const;
+  virtual Float_t  Unfold(Double_t eps, Int_t layer, Double_t *padSignal) const;
           Double_t GetCOG(Double_t signal[5]) const; 
   void             FillLUT();
           Double_t LUTposition(Int_t ilayer, Double_t ampL, Double_t ampC, Double_t ampR) const;
@@ -93,12 +107,9 @@ class AliTRDclusterizer : public TNamed
   UChar_t          GetPadStatus(UChar_t encoding) const;
   Int_t            GetCorruption(UChar_t encoding) const;
 
-  Bool_t           IsMaximum(const Int_t row, const Int_t col, const Int_t time, 
-                            UChar_t &pasStatus, Double_t *const Signals);
-  Bool_t           IsFivePadCluster(const Int_t row, const Int_t col, const Int_t time, 
-                                   Double_t *SignalsThisMax, Double_t *SignalsNeighbourMax, Double_t &ratio);
-  void             CreateCluster(const Int_t row, const Int_t col, const Int_t time, 
-                                const Double_t *const clusterSignal, const UChar_t padStatus);
+  Bool_t           IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Float_t *const Signals);       //for const correctness reasons not const parameters are given separately
+  Bool_t           IsFivePadCluster(const MaxStruct &ThisMax, const MaxStruct &NeighbourMax, Float_t &ratio); // ''
+  void             CreateCluster(const MaxStruct &Max); 
 
   const AliTRDReconstructor *fReconstructor;       //! reconstructor
   AliRunLoader        *fRunLoader;           //! Run Loader
@@ -139,20 +150,6 @@ class AliTRDclusterizer : public TNamed
   Int_t                fClusterROC;           // The index to the first cluster of a given ROC
   Int_t                firstClusterROC;       // The number of cluster in a given ROC
 
-  struct ClusterizerStruct
-  {
-    Int_t       Row;
-    Int_t       Col;
-    Int_t       Time;
-    UChar_t     padStatus;
-    Double_t    Signals[3];
-    ClusterizerStruct():Row(0),Col(0),Time(0),padStatus(0)
-      {}
-    ClusterizerStruct &operator=(const ClusterizerStruct &a)
-      {Row=a.Row; Col=a.Col; Time=a.Time; padStatus=a.padStatus;
-       memcpy(Signals, a.Signals, 3*sizeof(Double_t)); return *this;}
-  };
-
   ClassDef(AliTRDclusterizer,6)              //  TRD clusterfinder
 
 };
index f80989b..a6a41cc 100644 (file)
@@ -35,6 +35,7 @@
 #include "AliTRDdigit.h"
 #include "AliTRDgeometry.h"
 #include "AliTRDSignalIndex.h"
+#include "AliTRDReconstructor.h"
 
 ClassImp(AliTRDdigitsManager)
 
@@ -44,7 +45,7 @@ ClassImp(AliTRDdigitsManager)
   const Int_t AliTRDdigitsManager::fgkNDict = kNDict;
 
 //_____________________________________________________________________________
-AliTRDdigitsManager::AliTRDdigitsManager()
+AliTRDdigitsManager::AliTRDdigitsManager(const AliTRDReconstructor *const rec)
   :TObject()
   ,fEvent(0)
   ,fTree(0)
@@ -54,18 +55,23 @@ AliTRDdigitsManager::AliTRDdigitsManager()
   ,fUseDictionaries(kTRUE)
   ,fTreeD(0)
   ,fBranch(0)
+  ,fDets(AliTRDgeometry::Ndet())
+  ,fRawRec(kFALSE)
 {
   //
   // Default constructor
   //
+  
+  if(rec)
+    {
+      fDets=1;
+      fRawRec=kTRUE;
+    }
 
   for (Int_t iDict = 0; iDict < kNDict; iDict++) 
     {
       fDict[iDict] = NULL;
-    }
-  
-  fSignalIndexes = new TObjArray(AliTRDgeometry::Ndet());
-  
+    }  
 }
 
 //_____________________________________________________________________________
@@ -79,6 +85,8 @@ AliTRDdigitsManager::AliTRDdigitsManager(const AliTRDdigitsManager &m)
   ,fUseDictionaries(kTRUE)
   ,fTreeD(m.fTree)
   ,fBranch(m.fBranch)
+  ,fDets(m.fDets)
+  ,fRawRec(m.fRawRec)
 {
   //
   // AliTRDdigitsManager copy constructor
@@ -103,9 +111,12 @@ AliTRDdigitsManager::~AliTRDdigitsManager()
 
   for (Int_t iDict = 0; iDict < kNDict; iDict++) 
     {
-      fDict[iDict]->Delete();
-      delete fDict[iDict];
-      fDict[iDict] = NULL;
+      if(fDict[iDict])
+       {
+         fDict[iDict]->Delete();
+         delete fDict[iDict];
+         fDict[iDict] = NULL;
+       }
     }
 
   if (fSignalIndexes) 
@@ -149,6 +160,8 @@ void AliTRDdigitsManager::Copy(TObject &m) const
     }
   ((AliTRDdigitsManager &) m).fSignalIndexes   = fSignalIndexes;
   ((AliTRDdigitsManager &) m).fUseDictionaries = fUseDictionaries;
+  ((AliTRDdigitsManager &) m).fDets            = fDets;
+  ((AliTRDdigitsManager &) m).fRawRec           = fRawRec;
 
   TObject::Copy(m);
 
@@ -168,13 +181,9 @@ void AliTRDdigitsManager::CreateArrays()
          fDigits->Delete();                                
          delete fDigits;                                   
        }                                                    
-      fDigits = new TObjArray(AliTRDgeometry::Ndet());
-      for (Int_t index = 0; index < AliTRDgeometry::Ndet(); index++) 
-       {
-         AliTRDarraySignal *chamber= new AliTRDarraySignal();
-         chamber->SetNdet(index);
-         fDigits->AddAt(chamber,index);
-       } 
+      fDigits = new TObjArray(fDets);
+      for (Int_t index = 0; index < fDets; index++) 
+       fDigits->AddAt(new AliTRDarraySignal(),index);
     }
   else 
     {
@@ -183,13 +192,9 @@ void AliTRDdigitsManager::CreateArrays()
          fDigits->Delete();                                
          delete fDigits;                                   
        }                                                   
-      fDigits = new TObjArray(AliTRDgeometry::Ndet());    
-      for (Int_t index = 0; index < AliTRDgeometry::Ndet(); index++) 
-       {
-         AliTRDarrayADC *chamber= new AliTRDarrayADC();                                                 
-         chamber->SetNdet(index);
-         fDigits->AddAt(chamber,index);
-       }       
+      fDigits = new TObjArray(fDets);    
+      for (Int_t index = 0; index < fDets;index++) 
+       fDigits->AddAt(new AliTRDarrayADC(),index);
     }
 
   if (fUseDictionaries) 
@@ -201,25 +206,21 @@ void AliTRDdigitsManager::CreateArrays()
            delete fDict[iDict];                                    
          }
       for(Int_t iDict = 0; iDict < kNDict; iDict++)
-       {
-         fDict[iDict] = new TObjArray(AliTRDgeometry::Ndet());    
-       }
+       fDict[iDict] = new TObjArray(fDets);
 
       for (Int_t iDict = 0; iDict < kNDict; iDict++)
-       {
-         for (Int_t index = 0; index < AliTRDgeometry::Ndet(); index++) 
-           {
-             AliTRDarrayDictionary *dictio= new AliTRDarrayDictionary();
-             dictio->SetNdet(index);
-             fDict[iDict]->AddAt(dictio,index);
-           } 
-       }
+       for (Int_t index = 0; index < fDets; index++) 
+         fDict[iDict]->AddAt(new AliTRDarrayDictionary(),index);
     }
-
-  for (Int_t i = 0; i < AliTRDgeometry::Ndet(); i++) 
+  
+  if(fSignalIndexes)
     {
-      fSignalIndexes->AddLast(new AliTRDSignalIndex());
-    } 
+      fSignalIndexes->Delete();
+      delete fSignalIndexes;
+    }
+  fSignalIndexes = new TObjArray(fDets);
+  for (Int_t i = 0; i < fDets; i++) 
+    fSignalIndexes->AddLast(new AliTRDSignalIndex());
 
 }
 
@@ -230,65 +231,74 @@ void AliTRDdigitsManager::ResetArrays()
   // Reset the data arrays
   //
 
-  if (fDigits) 
+  if (fDigits)
     {
       fDigits->Delete();
       delete fDigits;
     }
-
-  for (Int_t iDict = 0; iDict < kNDict; iDict++)          
-    {                                                        
-      if (fDict[iDict])                                
-       {                                                    
-         fDict[iDict]->Delete();                      
-         delete fDict[iDict];                        
-       }
-    }
-
   if (fHasSDigits)
     {
-      fDigits = new TObjArray(AliTRDgeometry::Ndet());     
-      for (Int_t index = 0; index < AliTRDgeometry::Ndet(); index++) 
-       {
-         AliTRDarraySignal *chamber= new AliTRDarraySignal();
-         chamber->SetNdet(index);
-         fDigits->AddAt(chamber,index);
-       } 
+      fDigits = new TObjArray(fDets);     
+      for (Int_t index = 0; index < fDets; index++) 
+       fDigits->AddAt(new AliTRDarraySignal(),index);
     }
-  else 
+  else
+    {
+      fDigits = new TObjArray(fDets);
+      for (Int_t index = 0; index < fDets; index++)
+       fDigits->AddAt(new AliTRDarrayADC(),index);
+    }
+  
+  for (Int_t iDict = 0; iDict < kNDict; iDict++)
     {
-      fDigits = new TObjArray(AliTRDgeometry::Ndet());      
-      for (Int_t index = 0; index < AliTRDgeometry::Ndet(); index++) 
+      if (fDict[iDict])
        {
-         AliTRDarrayADC *chamber= new AliTRDarrayADC();
-         chamber->SetNdet(index);
-         fDigits->AddAt(chamber,index);
+         fDict[iDict]->Delete();
+         delete fDict[iDict];
+         fDict[iDict]=NULL;
        }
     }
-  
   if (fUseDictionaries) 
     {
       for(Int_t iDict = 0; iDict < kNDict; iDict++)
-       {
-         fDict[iDict] = new TObjArray(AliTRDgeometry::Ndet());  
-       }
+       fDict[iDict] = new TObjArray(fDets);
+      
       for (Int_t iDict = 0; iDict < kNDict; iDict++)
-       {
-         for (Int_t index = 0; index < AliTRDgeometry::Ndet(); index++) 
-           {
-             AliTRDarrayDictionary *dictio= new AliTRDarrayDictionary();   
-             dictio->SetNdet(index);                                       
-             fDict[iDict]->AddAt(dictio,index);                            
-           }
-       }
+       for (Int_t index = 0; index < fDets; index++)
+         fDict[iDict]->AddAt(new AliTRDarrayDictionary(),index);
     }
-
-  for (Int_t i = 0; i < AliTRDgeometry::Ndet(); i++) 
+  
+  if(fSignalIndexes)
     {
-      AliTRDSignalIndex *idx = (AliTRDSignalIndex *) fSignalIndexes->At(i);
-      if (idx) idx->Reset();
+      fSignalIndexes->Delete();
+      delete fSignalIndexes;
     }
+  fSignalIndexes = new TObjArray(fDets);
+  for (Int_t i = 0; i < fDets; i++)
+    fSignalIndexes->AddLast(new AliTRDSignalIndex());
+}
+
+//_____________________________________________________________________________
+void AliTRDdigitsManager::ResetArrays(Int_t det)
+{
+  Int_t RecoDet = fRawRec ? 0 : det;
+
+  RemoveDigits(RecoDet);
+  RemoveDictionaries(RecoDet);
+  RemoveIndexes(RecoDet);
+
+  if (fHasSDigits)
+    fDigits->AddAt(new AliTRDarraySignal(),RecoDet);
+  else
+    fDigits->AddAt(new AliTRDarrayADC(),RecoDet);
 
+  if (fUseDictionaries) 
+    {
+      for (Int_t iDict = 0; iDict < kNDict; iDict++)
+       fDict[iDict]->AddAt(new AliTRDarrayDictionary(),RecoDet);
+    }
+  
+  fSignalIndexes->AddAt(new AliTRDSignalIndex(),RecoDet);
 }
 
 //_____________________________________________________________________________
@@ -298,10 +308,7 @@ Short_t AliTRDdigitsManager::GetDigitAmp(Int_t row, Int_t col,Int_t time, Int_t
   // Returns the amplitude of a digit
   //
 
-  if (!GetDigits(det)) 
-    {
-      return 0;
-    }
+  if (!GetDigits(det)) return 0;
   
   return ((Short_t) ((AliTRDarrayADC *) GetDigits(det))->GetDataB(row,col,time));
 
@@ -314,10 +321,7 @@ UChar_t AliTRDdigitsManager::GetPadStatus(Int_t row, Int_t col, Int_t time, Int_
   // Returns the pad status for the requested pad
   //
        
-  if (!GetDigits(det)) 
-    {
-      return 0;
-    }
+  if (!GetDigits(det)) return 0;
 
   return ((UChar_t) ((AliTRDarrayADC *) GetDigits(det))->GetPadStatus(row,col,time));
  
@@ -551,6 +555,8 @@ AliTRDarrayADC *AliTRDdigitsManager::GetDigits(Int_t det) const
   // Returns the digits array for one detector
   //
 
+  Int_t RecoDet = fRawRec ? 0 : det;
+
   if (!fDigits)   
     {
       return 0x0;
@@ -558,7 +564,8 @@ AliTRDarrayADC *AliTRDdigitsManager::GetDigits(Int_t det) const
 
   if (!fHasSDigits)
     {
-      return (AliTRDarrayADC *) fDigits->At(det); 
+      ((AliTRDarrayADC *) fDigits->At(RecoDet))->SetNdet(det);
+      return (AliTRDarrayADC *) fDigits->At(RecoDet); 
     }
   else
     {
@@ -575,6 +582,8 @@ AliTRDarraySignal *AliTRDdigitsManager::GetSDigits(Int_t det) const
   // Returns the sdigits array for one detector
   //
 
+  Int_t RecoDet = fRawRec ? 0 : det;
+
   if (!fDigits)   
     {
       //      AliDebug(1,"NO FDIGITS!");       
@@ -583,7 +592,8 @@ AliTRDarraySignal *AliTRDdigitsManager::GetSDigits(Int_t det) const
 
   if (fHasSDigits)
     {
-      return (AliTRDarraySignal *) fDigits->At(det);
+      ((AliTRDarraySignal *) fDigits->At(RecoDet))->SetNdet(det);
+      return (AliTRDarraySignal *) fDigits->At(RecoDet);
     }
   else
     {
@@ -601,12 +611,15 @@ AliTRDarrayDictionary *AliTRDdigitsManager::GetDictionary(Int_t det
   // Returns the dictionary for one detector
   //
 
+  Int_t RecoDet = fRawRec ? 0 : det;
+
   if (fUseDictionaries == kFALSE)
     {
       return 0x0;
     }
 
-  return (AliTRDarrayDictionary *) fDict[i]->At(det);
+  ((AliTRDarrayDictionary *) fDigits->At(RecoDet))->SetNdet(det);
+  return (AliTRDarrayDictionary *) fDict[i]->At(RecoDet);
   
 }
 
@@ -633,7 +646,9 @@ AliTRDSignalIndex *AliTRDdigitsManager::GetIndexes(Int_t det)
   // Returns indexes of active pads
   //
 
-  return (AliTRDSignalIndex *) fSignalIndexes->At(det);
+  Int_t RecoDet = fRawRec ? 0 : det;
+
+  return (AliTRDSignalIndex *) fSignalIndexes->At(RecoDet);
 
 }
 
@@ -644,16 +659,18 @@ void AliTRDdigitsManager::RemoveDigits(Int_t det)
    // Clear memory at det for Digits
    //
 
-  if (fDigits->At(det))
+  Int_t RecoDet = fRawRec ? 0 : det;
+
+  if (fDigits->At(RecoDet))
     {
       if (fHasSDigits) 
         {
-          AliTRDarraySignal *arr = (AliTRDarraySignal *) fDigits->RemoveAt(det);
+          AliTRDarraySignal *arr = (AliTRDarraySignal *) fDigits->RemoveAt(RecoDet);
           delete arr;
        }
       else 
         {
-          AliTRDarrayADC    *arr = (AliTRDarrayADC *)    fDigits->RemoveAt(det);
+          AliTRDarrayADC    *arr = (AliTRDarrayADC *)    fDigits->RemoveAt(RecoDet);
           delete arr;
        }
     }
@@ -667,6 +684,8 @@ void AliTRDdigitsManager::RemoveDictionaries(Int_t det)
   // Clear memory
   //
 
+  Int_t RecoDet = fRawRec ? 0 : det;
+
   if (fUseDictionaries == kFALSE) 
     {
       return;
@@ -674,23 +693,43 @@ void AliTRDdigitsManager::RemoveDictionaries(Int_t det)
 
   for (Int_t i = 0; i < kNDict; i++) 
     {
-      if (fDict[i]->At(det))
+      if (fDict[i]->At(RecoDet))
        {
-         AliTRDarrayDictionary *arr = (AliTRDarrayDictionary *) fDict[i]->RemoveAt(det);
+         AliTRDarrayDictionary *arr = (AliTRDarrayDictionary *) fDict[i]->RemoveAt(RecoDet);
           delete arr;
        }
     }
 
 }
 
+//_____________________________________________________________________________
+void AliTRDdigitsManager::RemoveIndexes(Int_t det) 
+{
+   // 
+   // Clear memory
+   //
+
+  Int_t RecoDet = fRawRec ? 0 : det;
+
+  if (fSignalIndexes->At(RecoDet))
+    {
+      AliTRDSignalIndex *arr = (AliTRDSignalIndex *) fSignalIndexes->RemoveAt(RecoDet);
+      delete arr;
+    }
+
+}
+
+
 //_____________________________________________________________________________
 void AliTRDdigitsManager::ClearIndexes(Int_t det) 
 {
   // 
   // Clear memory
   //
+  
+  Int_t RecoDet = fRawRec ? 0 : det;
 
-  ((AliTRDSignalIndex *) fSignalIndexes->At(det))->ClearAll();  
+  ((AliTRDSignalIndex *) fSignalIndexes->At(RecoDet))->ClearAll();  
 
 }
 
@@ -790,7 +829,7 @@ Bool_t AliTRDdigitsManager::LoadArray(TObjArray *object
 
   // Loop through all detectors and read them from the tree
   Bool_t status = kTRUE;
-  for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) 
+  for (Int_t iDet = 0; iDet < fDets; iDet++) 
     {
       if(fHasSDigits)
        {
@@ -850,7 +889,7 @@ Bool_t AliTRDdigitsManager::LoadArrayDict(TObjArray *object
 
   // Loop through all detectors and read them from the tree
   Bool_t status = kTRUE;
-  for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) 
+  for (Int_t iDet = 0; iDet < fDets; iDet++) 
     {
       AliTRDarrayDictionary *dataArray = (AliTRDarrayDictionary *) object->At(iDet);
       if (!dataArray) 
@@ -895,7 +934,7 @@ Bool_t AliTRDdigitsManager::StoreArray(TObjArray *array1
 
   // Loop through all detectors and fill them into the tree
   Bool_t status = kTRUE;
-  for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) 
+  for (Int_t iDet = 0; iDet < fDets; iDet++) 
     {
       if (fHasSDigits)
        {
@@ -955,7 +994,7 @@ Bool_t AliTRDdigitsManager::StoreArrayDict(TObjArray *array3
 
   // Loop through all detectors and fill them into the tree
   Bool_t status = kTRUE;
-  for (Int_t iDet = 0; iDet < AliTRDgeometry::Ndet(); iDet++) 
+  for (Int_t iDet = 0; iDet < fDets; iDet++) 
     {
       const AliTRDarrayDictionary *kDataArray = (AliTRDarrayDictionary *) array3->At(iDet);
       if (!kDataArray) 
index 10db6a9..33c01e1 100644 (file)
@@ -20,6 +20,7 @@ class AliTRDSignalIndex;
 class AliTRDarrayADC;  
 class AliTRDarraySignal; 
 class AliTRDarrayDictionary;
+class AliTRDReconstructor;
 
 class AliTRDdigitsManager : public TObject {
 
@@ -27,7 +28,7 @@ class AliTRDdigitsManager : public TObject {
 
   enum { kNDict = 3 };
 
-  AliTRDdigitsManager();
+  AliTRDdigitsManager(const AliTRDReconstructor *const rec = 0x0);  //if rec is given, we are reading raw data, so the TObjArrays may (and will) contain only one entry
   AliTRDdigitsManager(const AliTRDdigitsManager &m);
   virtual ~AliTRDdigitsManager();
   AliTRDdigitsManager &operator=(const AliTRDdigitsManager &m);
@@ -36,6 +37,7 @@ class AliTRDdigitsManager : public TObject {
 
   virtual void                CreateArrays();
   virtual void                ResetArrays();
+  virtual void                ResetArrays(Int_t det);
   virtual Bool_t              BuildIndexes(Int_t det);
 
   virtual Bool_t              MakeBranch(TTree *tree);
@@ -65,6 +67,7 @@ class AliTRDdigitsManager : public TObject {
 
   void                        RemoveDigits(Int_t det);
   void                        RemoveDictionaries(Int_t det);
+  void                        RemoveIndexes(Int_t det);
   void                        ClearIndexes(Int_t det);
   
   Int_t                       GetTrack(Int_t track, AliTRDdigit *digit) const;
@@ -88,6 +91,9 @@ class AliTRDdigitsManager : public TObject {
   Bool_t              fUseDictionaries;    //  Use dictionaries or not (case of real data)
   TTree              *fTreeD;              //  Tree with detector objects
   TBranch            *fBranch;             //  Branchaddress
+  Int_t               fDets;               //  No of Detectors
+  Bool_t              fRawRec;             //  Reconstruct from raw data. If its kTRUE then the TObjArrays have only one entry.
+                                           //  If kFALSE then they have (AliTRDgeometry::Ndet()) entries (default).
 
   ClassDef(AliTRDdigitsManager,7)          //  Manages the TRD digits