]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDclusterizer.cxx
Fix in procedure for vertex recalculation
[u/mrichter/AliRoot.git] / TRD / AliTRDclusterizer.cxx
index ab89b678bf18c9ef57ebc74a2e46ad122e02159f..60bad65c4d9d82b9b0264682e782f03a56491e3a 100644 (file)
@@ -41,7 +41,6 @@
 #include "AliTRDcalibDB.h"
 #include "AliTRDtransform.h"
 #include "AliTRDSignalIndex.h"
-#include "AliTRDrawStreamBase.h"
 #include "AliTRDrawStream.h"
 #include "AliTRDfeeParam.h"
 #include "AliTRDtrackletWord.h"
@@ -62,6 +61,7 @@ AliTRDclusterizer::AliTRDclusterizer(const AliTRDReconstructor *const rec)
   ,fClusterTree(NULL)
   ,fRecPoints(NULL)
   ,fTracklets(NULL)
+  ,fTracks(NULL)
   ,fTrackletTree(NULL)
   ,fDigitsManager(new AliTRDdigitsManager())
   ,fTrackletContainer(NULL)
@@ -120,6 +120,7 @@ AliTRDclusterizer::AliTRDclusterizer(const Text_t *name
   ,fClusterTree(NULL)
   ,fRecPoints(NULL)
   ,fTracklets(NULL)
+  ,fTracks(NULL)
   ,fTrackletTree(NULL)
   ,fDigitsManager(new AliTRDdigitsManager())
   ,fTrackletContainer(NULL)
@@ -171,6 +172,7 @@ AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c)
   ,fClusterTree(NULL)
   ,fRecPoints(NULL)
   ,fTracklets(NULL)
+  ,fTracks(NULL)
   ,fTrackletTree(NULL)
   ,fDigitsManager(NULL)
   ,fTrackletContainer(NULL)
@@ -227,18 +229,16 @@ AliTRDclusterizer::~AliTRDclusterizer()
     delete fTracklets;
   }
 
+  if (fTracks){
+    fTracks->Delete();
+    delete fTracks;
+  }
+
   if (fDigitsManager) {
     delete fDigitsManager;
     fDigitsManager = NULL;
   }
 
-  if (fTrackletContainer){
-    delete [] fTrackletContainer[0];
-    delete [] fTrackletContainer[1];
-    delete [] fTrackletContainer;
-    fTrackletContainer = NULL;
-  }
-
   if (fTransform){
     delete fTransform;
     fTransform = NULL;
@@ -278,7 +278,6 @@ void AliTRDclusterizer::Copy(TObject &c) const
   ((AliTRDclusterizer &) c).fRecPoints     = NULL;  
   ((AliTRDclusterizer &) c).fTrackletTree  = NULL;
   ((AliTRDclusterizer &) c).fDigitsManager = NULL;
-  ((AliTRDclusterizer &) c).fTrackletContainer = NULL;
   ((AliTRDclusterizer &) c).fRawVersion    = fRawVersion;
   ((AliTRDclusterizer &) c).fTransform     = NULL;
   ((AliTRDclusterizer &) c).fDigits      = NULL;
@@ -436,48 +435,6 @@ Bool_t AliTRDclusterizer::WriteClusters(Int_t det)
   return kTRUE;  
 }
 
-//_____________________________________________________________________________
-Bool_t AliTRDclusterizer::WriteTracklets(Int_t det)
-{
-  //
-  // Write the raw data tracklets into seperate file
-  //
-
-  UInt_t **leaves = new UInt_t *[2];
-  for (Int_t i=0; i<2 ;i++){
-    leaves[i] = new UInt_t[258];
-    leaves[i][0] = det; // det
-    leaves[i][1] = i;   // side
-    memcpy(leaves[i]+2, fTrackletContainer[i], sizeof(UInt_t) * 256);
-  }
-
-  if (!fTrackletTree){
-    AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
-    dl->MakeTree();
-    fTrackletTree = dl->Tree();
-  }
-
-  TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
-  if (!trkbranch) {
-    trkbranch = fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
-  }
-
-  for (Int_t i=0; i<2; i++){
-    if (leaves[i][2]>0) {
-      trkbranch->SetAddress(leaves[i]);
-      fTrackletTree->Fill();
-    }
-  }
-
-//  AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets"); //jkl: wrong here
-//  dl->WriteData("OVERWRITE"); //jkl: wrong here
-  //dl->Unload();
-  delete [] leaves;
-
-  return kTRUE;
-
-}
-
 //_____________________________________________________________________________
 Bool_t AliTRDclusterizer::ReadDigits()
 {
@@ -625,44 +582,28 @@ Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
       trklTreeLoader->MakeTree();
   }
 
-  // tracklet container for raw tracklet writing
-  if (!fTrackletContainer && ( fReconstructor->IsWritingTracklets() || fReconstructor->IsProcessingTracklets() )) {
-    // maximum tracklets for one HC
-    const Int_t kTrackletChmb=256;
-    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)
-    fRawStream = AliTRDrawStreamBase::GetRawStream(rawReader);
+    fRawStream = new AliTRDrawStream(rawReader);
   else
     fRawStream->SetReader(rawReader);
 
-  if(fReconstructor->IsHLT()){
-    if(fRawStream->InheritsFrom(AliTRDrawStream::Class()))
-      ((AliTRDrawStream*)fRawStream)->DisableErrorStorage();
-    else{
-      fRawStream->SetSharedPadReadout(kFALSE);
-      fRawStream->SetNoErrorWarning();
-    }
-  }
+  //if(fReconstructor->IsHLT()){
+    fRawStream->DisableErrorStorage();
+  //}
+
+  // register tracklet array for output
+  if (fReconstructor->IsProcessingTracklets())
+    fRawStream->SetTrackletArray(TrackletsArray());
+  fRawStream->SetTrackArray(TracksArray());
 
-  AliDebug(1,Form("Stream version: %s", fRawStream->IsA()->GetName()));
-  
   UInt_t det = 0;
-  while ((det = fRawStream->NextChamber(fDigitsManager,fTrackletContainer)) < AliTRDgeometry::kNdet){
+  while ((det = fRawStream->NextChamber(fDigitsManager)) < AliTRDgeometry::kNdet){
     if (fDigitsManager->GetIndexes(det)->HasEntry()){
       MakeClusters(det);
       fDigitsManager->ClearArrays(det);
     }
-    
-    if (!fReconstructor->IsWritingTracklets()) continue;
-    if (*(fTrackletContainer[0]) > 0 || *(fTrackletContainer[1]) > 0) WriteTracklets(det);
   }
-  
+
   if (fReconstructor->IsWritingTracklets()) {
     if (AliDataLoader *trklLoader = AliRunLoader::Instance()->GetLoader("TRDLoader")->GetDataLoader("tracklets")) {
       if (trklLoader) {
@@ -673,13 +614,6 @@ Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
     }
   }
 
-  if (fTrackletContainer){
-    delete [] fTrackletContainer[0];
-    delete [] fTrackletContainer[1];
-    delete [] fTrackletContainer;
-    fTrackletContainer = NULL;
-  }
-
   if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
 
   if(!TestBit(knewDM)){
@@ -801,9 +735,6 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
   Int_t    iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
   fVolid      = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule); 
 
-  if(fReconstructor->IsProcessingTracklets() && fTrackletContainer)
-    AddTrackletsToArray();
-
   fColMax    = fDigits->GetNcol();
   fTimeTotal = fDigitsManager->GetDigitsParam()->GetNTimeBins(det);
 
@@ -855,7 +786,9 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
   SetBit(kGAUS, recoParam->UseGAUS());
 
   // Apply the gain and the tail cancelation via digital filter
-  if(recoParam->UseTailCancelation()) TailCancelation(recoParam);
+  // Use the configuration from the DCS to find out whether online 
+  // tail cancellation was applied
+  if(!calibration->HasOnlineTailCancellation()) TailCancelation(recoParam);
 
   MaxStruct curr, last;
   Int_t nMaximas = 0, nCorrupted = 0;
@@ -884,7 +817,7 @@ Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
       << "NCorrupted=" << nCorrupted
       << "\n";
   }
-  if (TestBit(kLabels)) AddLabels();
+  // if (TestBit(kLabels)) AddLabels();
 
   return kTRUE;
 
@@ -898,8 +831,10 @@ Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Sh
   // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
   //
 
+  Float_t tmp(0.), kMaxShortVal(32767.); // protect against data overflow due to wrong gain calibration
   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);
+  tmp = (fDigits->GetData(Max.row, Max.col, Max.time) - fBaseline) / gain + 0.5f;
+  Signals[1] = (Short_t)TMath::Min(tmp, kMaxShortVal);
   if(Signals[1] <= fMaxThresh) return kFALSE;
 
   if(Max.col < 1 || Max.col + 1 >= fColMax) return kFALSE;
@@ -913,10 +848,17 @@ 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);
+  Short_t signal(0);
+  if((signal = fDigits->GetData(Max.row, Max.col-1, Max.time))){
+    gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
+    tmp = (signal - fBaseline) / gain + 0.5f;
+    Signals[0] = (Short_t)TMath::Min(tmp, kMaxShortVal);
+  } else Signals[0] = 0;
+  if((signal = fDigits->GetData(Max.row, Max.col+1, Max.time))){
+    gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+1,Max.row);
+    tmp = (signal - fBaseline) / gain + 0.5f;
+    Signals[2] = (Short_t)TMath::Min(tmp, kMaxShortVal);
+  } else Signals[2] = 0;
 
   if(!(status[0] | status[1] | status[2])) {//all pads are good
     if ((Signals[2] <= Signals[1]) && (Signals[0] <  Signals[1])) {
@@ -1038,42 +980,41 @@ void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
 //_____________________________________________________________________________
 void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
 {
+// Calculate number of pads/cluster and
+// ADC signals at position 0, 1, 5 and 6
+
+  Float_t tmp(0.), kMaxShortVal(32767.); // protect against data overflow due to wrong gain calibration
+  Float_t gain(1.); Short_t signal(0);
+  // Store the amplitudes of the pads in the cluster for later analysis
+  // and check whether one of these pads is masked in the database
+  signals[3]=Max.signals[1];
+  Int_t ipad(1), jpad(0);
   // Look to the right
-  Int_t ii = 1;
-  while (fDigits->GetData(Max.row, Max.col-ii, Max.time) >= fSigThresh) {
+  while((jpad = Max.col-ipad)){
+    if(!(signal = fDigits->GetData(Max.row, jpad, Max.time))) break; // empty digit !
+    gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row);
+    tmp = (signal - fBaseline) / gain + 0.5f;
+    signal = (Short_t)TMath::Min(tmp, kMaxShortVal);
+    if(signal<fSigThresh) break; // signal under threshold
     nPadCount++;
-    ii++;
-    if (Max.col < ii) break;
+    if(ipad<=3) signals[3 - ipad] = signal;
+    ipad++;
   }
+  ipad=1;
   // Look to the left
-  ii = 1;
-  while (fDigits->GetData(Max.row, Max.col+ii, Max.time) >= fSigThresh) {
+  while((jpad = Max.col+ipad)<fColMax){
+    if(!(signal = fDigits->GetData(Max.row, jpad, Max.time))) break; // empty digit !
+    gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row);
+    tmp = (signal - fBaseline) / gain + 0.5f;
+    signal = (Short_t)TMath::Min(tmp, kMaxShortVal);
+    if(signal<fSigThresh) break; // signal under threshold
     nPadCount++;
-    ii++;
-    if (Max.col+ii >= fColMax) break;
+    if(ipad<=3) signals[3 + ipad] = signal;
+    ipad++;
   }
 
-  // Store the amplitudes of the pads in the cluster for later analysis
-  // and check whether one of these pads is masked in the database
-  signals[2]=Max.signals[0];
-  signals[3]=Max.signals[1];
-  signals[4]=Max.signals[2];
-  Float_t gain;
-  for(Int_t i = 0; i<2; i++)
-    {
-      if(Max.col+i >= 3){
-       gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-3+i,Max.row);
-       signals[i] = (Short_t)((fDigits->GetData(Max.row, Max.col-3+i, Max.time) - fBaseline) / gain + 0.5f);
-      }
-      if(Max.col+3-i < fColMax){
-       gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+3-i,Max.row);
-       signals[6-i] = (Short_t)((fDigits->GetData(Max.row, Max.col+3-i, Max.time) - fBaseline) / gain + 0.5f);
-      }
-    }
-  /*for (Int_t jPad = Max.Col-3; jPad <= Max.Col+3; jPad++) {
-    if ((jPad >= 0) && (jPad < fColMax))
-      signals[jPad-Max.Col+3] = TMath::Nint(fDigits->GetData(Max.Row,jPad,Max.Time));
-      }*/
+  AliDebug(4, Form("Signals[%3d %3d %3d %3d %3d %3d %3d] Npads[%d]."
+    , signals[0], signals[1], signals[2], signals[3], signals[4], signals[5], signals[6], nPadCount));
 }
 
 //_____________________________________________________________________________
@@ -1088,27 +1029,6 @@ void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
   new((*RecPoints())[n]) AliTRDcluster(*cluster);
 }
 
-//_____________________________________________________________________________
-void AliTRDclusterizer::AddTrackletsToArray()
-{
-  //
-  // Add the online tracklets of this chamber to the array
-  //
-
-  UInt_t* trackletword;
-  for(Int_t side=0; side<2; side++)
-    {
-      Int_t trkl=0;
-      trackletword=fTrackletContainer[side];
-      while(trackletword[trkl]>0){
-       Int_t n = TrackletsArray()->GetEntriesFast();
-       AliTRDtrackletWord tmp(trackletword[trkl]);
-       new((*TrackletsArray())[n]) AliTRDcluster(&tmp,fDet,fVolid);
-       trkl++;
-      }
-    }
-}
-
 //_____________________________________________________________________________
 Bool_t AliTRDclusterizer::AddLabels()
 {
@@ -1242,13 +1162,15 @@ void AliTRDclusterizer::TailCancelation(const AliTRDrecoParam* const recoParam)
   // Applies the tail cancelation
   //
 
+  Int_t nexp = recoParam->GetTCnexp();
+  if(!nexp) return;
+  
   Int_t iRow  = 0;
   Int_t iCol  = 0;
   Int_t iTime = 0;
 
   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))
     {
       // if corrupted then don't make the tail cancallation
@@ -1263,14 +1185,78 @@ void AliTRDclusterizer::TailCancelation(const AliTRDrecoParam* const recoParam)
       }
       
       // Apply the tail cancelation via the digital filter
-      DeConvExp(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp);
-      
+      //DeConvExp(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp);
+      ApplyTCTM(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp);
     } // while irow icol
 
   return;
 
 }
 
+
+//_____________________________________________________________________________
+void AliTRDclusterizer::ApplyTCTM(Short_t *const arr, const Int_t nTime, const Int_t nexp) 
+{
+  //
+  // Steer tail cancellation
+  //
+
+
+  switch(nexp) {
+  case 1:
+  case 2:
+    DeConvExp(arr,nTime,nexp);
+    break;
+  case -1:
+    ConvExp(arr,nTime);
+    break;
+  case -2:
+    DeConvExp(arr,nTime,1);
+    ConvExp(arr,nTime);
+    break;
+  default:
+    break;
+  }
+}
+
+
+//_____________________________________________________________________________
+void AliTRDclusterizer::ConvExp(Short_t *const arr, const Int_t nTime)
+{
+  //
+  // Tail maker
+  //
+
+  // Initialization (coefficient = alpha, rates = lambda)
+  Float_t slope = 1.0;
+  Float_t coeff = 0.5;
+  Float_t rate;
+
+  Double_t par[4];
+  fReconstructor->GetRecoParam()->GetTCParams(par);
+  slope = par[1];
+  coeff = par[3];  
+
+  Double_t dt = 0.1;
+
+  rate = TMath::Exp(-dt/(slope));
+   
+  Float_t reminder =  .0;
+  Float_t correction = 0.0;
+  Float_t result     = 0.0;
+
+  for (int i = nTime-1; i >= 0; i--) {
+
+    result = arr[i] + correction - fBaseline;    // No rescaling
+    arr[i] = (Short_t)(result + fBaseline + 0.5f);
+
+    correction = 0.0;
+    
+    correction += reminder = rate * (reminder + coeff * result);
+  }
+}
+
+
 //_____________________________________________________________________________
 void AliTRDclusterizer::DeConvExp(Short_t *const arr, const Int_t nTime, const Int_t nexp)
 {
@@ -1367,11 +1353,11 @@ TClonesArray *AliTRDclusterizer::RecPoints()
 TClonesArray *AliTRDclusterizer::TrackletsArray() 
 {
   //
-  // Returns the list of rec points
+  // Returns the array of on-line tracklets
   //
 
   if (!fTracklets && fReconstructor->IsProcessingTracklets()) {
-    fTracklets = new TClonesArray("AliTRDcluster", 2*MAXTRACKLETSPERHC);
+    fTracklets = new TClonesArray("AliTRDtrackletWord", 2*MAXTRACKLETSPERHC);
     //SetClustersOwner(kTRUE);
     //AliTRDReconstructor::SetTracklets(0x0);
   }
@@ -1379,3 +1365,13 @@ TClonesArray *AliTRDclusterizer::TrackletsArray()
 
 }
 
+//_____________________________________________________________________________
+TClonesArray* AliTRDclusterizer::TracksArray()
+{
+  // return array of GTU tracks (create TClonesArray if necessary)
+
+  if (!fTracks) {
+    fTracks = new TClonesArray("AliESDTrdTrack",100);
+  }
+  return fTracks;
+}