]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCclustererMI.cxx
Fixes for CAF, etc
[u/mrichter/AliRoot.git] / TPC / AliTPCclustererMI.cxx
index 95170c532976d37b2a2b8c0fe7b6f5f5abdeca3d..8b481f7b859fcd6a1b14c4fe3753599f10582d6d 100644 (file)
 //  1. The Input data for reconstruction - Options
 //      1.a Simulated data  - TTree - invoked Digits2Clusters()
 //      1.b Raw data        - Digits2Clusters(AliRawReader* rawReader); 
+//      1.c HLT clusters    - Digits2Clusters and Digits2Clusters(AliRawReader* rawReader)
+//                            invoke ReadHLTClusters()
+//
+//      fUseHLTClusters     - switches between different inputs
+//                            1 -> only TPC raw/sim data
+//                            2 -> if present TPC raw/sim data, otherwise HLT clusters
+//                            3 -> only HLT clusters
+//                            4 -> if present HLT clusters, otherwise TPC raw/sim data
 //
 //  2. The Output data
 //      2.a TTree with clusters - if  SetOutput(TTree * tree) invoked
@@ -48,6 +56,8 @@
 #include <TRandom.h>
 #include <TTree.h>
 #include <TTreeStream.h>
+#include "TSystem.h"
+#include "TClass.h"
 
 #include "AliDigits.h"
 #include "AliLoader.h"
@@ -72,6 +82,8 @@
 #include "AliTPCTransform.h"
 #include "AliTPCclustererMI.h"
 
+using std::cerr;
+using std::endl;
 ClassImp(AliTPCclustererMI)
 
 
@@ -108,9 +120,11 @@ AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoPar
   fRecoParam(0),
   fBDumpSignal(kFALSE),
   fBClonesArray(kFALSE),
+  fUseHLTClusters(4),
   fAllBins(NULL),
   fAllSigBins(NULL),
-  fAllNSigBins(NULL)
+  fAllNSigBins(NULL),
+  fHLTClusterAccess(NULL)
 {
   //
   // COSNTRUCTOR
@@ -154,58 +168,7 @@ AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoPar
     fAllNSigBins[iRow]=0;
   }
 }
-//______________________________________________________________
-AliTPCclustererMI::AliTPCclustererMI(const AliTPCclustererMI &param)
-              :TObject(param),
-  fBins(0),
-  fSigBins(0),
-  fNSigBins(0),
-  fLoop(0),
-  fMaxBin(0),
-  fMaxTime(0),
-  fMaxPad(0),
-  fSector(-1),
-  fRow(-1),
-  fSign(0),
-  fRx(0),
-  fPadWidth(0),
-  fPadLength(0),
-  fZWidth(0),
-  fPedSubtraction(kFALSE),
-  fEventHeader(0),
-  fTimeStamp(0),
-  fEventType(0),
-  fInput(0),
-  fOutput(0),
-  fOutputArray(0),
-  fOutputClonesArray(0),
-  fRowCl(0),
-  fRowDig(0),
-  fParam(0),
-  fNcluster(0),
-  fNclusters(0),
-  fDebugStreamer(0),
-  fRecoParam(0),
-  fBDumpSignal(kFALSE),
-  fBClonesArray(kFALSE),
-  fAllBins(NULL),
-  fAllSigBins(NULL),
-  fAllNSigBins(NULL)
-{
-  //
-  // dummy
-  //
-  fMaxBin = param.fMaxBin;
-}
-//______________________________________________________________
-AliTPCclustererMI & AliTPCclustererMI::operator =(const AliTPCclustererMI & param)
-{
-  //
-  // assignment operator - dummy
-  //
-  fMaxBin=param.fMaxBin;
-  return (*this);
-}
+
 //______________________________________________________________
 AliTPCclustererMI::~AliTPCclustererMI(){
   //
@@ -235,6 +198,7 @@ AliTPCclustererMI::~AliTPCclustererMI(){
   delete [] fAllBins;
   delete [] fAllSigBins;
   delete [] fAllNSigBins;
+  if (fHLTClusterAccess) delete fHLTClusterAccess;
 }
 
 void AliTPCclustererMI::SetInput(TTree * tree)
@@ -307,11 +271,25 @@ void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Float_t *bins, UInt_t /*m*
 AliTPCclusterMI &c) 
 {
   //
+  //  Make cluster: characterized by position ( mean-  COG) , shape (RMS) a charge, QMax and Q tot
+  //  Additional correction:
+  //  a) To correct for charge below threshold, in the +1 neghborhood to the max charge charge 
+  //       is extrapolated using gaussian approximation assuming given cluster width.. 
+  //       Additional empirical factor is used to account for the charge fluctuation (kVirtualChargeFactor). 
+  //       Actual value of the  kVirtualChargeFactor should obtained minimimizing residuals between the cluster
+  //       and track interpolation.
+  //  b.) For space points with extended shape (in comparison with expected using parameterization) clusters are 
+  //      unfoded     
+  //  
+  //  NOTE. Actual/Empirical  values for correction are hardwired in the code.
+  //
+  // Input paramters for function:
   //  k    - Make cluster at position k  
   //  bins - 2 D array of signals mapped to 1 dimensional array - 
   //  max  - the number of time bins er one dimension
-  //  c    - refernce to cluster to be filled
+  //  c    - reference to cluster to be filled
   //
+  Double_t kVirtualChargeFactor=0.5;
   Int_t i0=k/max;  //central pad
   Int_t j0=k%max;  //central time bin
 
@@ -337,7 +315,7 @@ AliTPCclusterMI &c)
        Float_t ratio = TMath::Exp(-1.2*TMath::Abs(di)/sigmay2)*TMath::Exp(-1.2*TMath::Abs(dj)/sigmaz2);
        amp = ((matrix[2][0]-2)*(matrix[2][0]-2)/(matrix[-di+2][-dj]+2))*ratio;
        if (amp>2) amp = 2;
-       vmatrix[2+di][2+dj]=amp;
+       vmatrix[2+di][2+dj]= kVirtualChargeFactor*amp;
        vmatrix[2+2*di][2+2*dj]=0;
        if ( (di*dj)!=0){       
          //DIAGONAL ELEMENTS
@@ -647,7 +625,6 @@ void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * /*matrix*/, Int
     c.SetType(-(c.GetType()+3));  //edge clusters
   }
   if (fLoop==2) c.SetType(100);
-  if (!AcceptCluster(&c)) return;
 
   // select output 
   TClonesArray * arr = 0;
@@ -709,7 +686,40 @@ void AliTPCclustererMI::Digits2Clusters()
     fParam->Dump();
     fRecoParam->Dump();
   }
+  fRowDig = NULL;
+
+  //-----------------------------------------------------------------
+  // Use HLT clusters
+  //-----------------------------------------------------------------
+  if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
+    AliInfo("Using HLT clusters for TPC off-line reconstruction");
+    fZWidth = fParam->GetZWidth();
+    Int_t iResult = ReadHLTClusters();
+
+    // HLT clusters present
+    if (iResult >= 0 && fNclusters > 0)
+      return; 
 
+    // HLT clusters not present
+    if (iResult < 0 || fNclusters == 0) {
+      if (fUseHLTClusters == 3) { 
+       AliError("No HLT clusters present, but requested.");
+       return;
+      }
+      else {
+       AliInfo("Now trying to read from TPC RAW");
+      }
+    }
+    // Some other problem during cluster reading
+    else {
+      AliWarning("Some problem while unpacking of HLT clusters.");
+      return;
+    }
+  } // if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
+
+  //-----------------------------------------------------------------
+  // Run TPC off-line clusterer
+  //-----------------------------------------------------------------
   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
   AliSimDigits digarr, *dummy=&digarr;
@@ -741,7 +751,7 @@ void AliTPCclustererMI::Digits2Clusters()
     fZWidth = fParam->GetZWidth();
     if (fSector < kNIS) {
       fMaxPad = fParam->GetNPadsLow(row);
-      fSign = (fSector < kNIS/2) ? 1 : -1;
+      fSign =  (fSector < kNIS/2) ? 1 : -1;
       fPadLength = fParam->GetPadPitchLength(fSector,row);
       fPadWidth = fParam->GetPadPitchWidth();
     } else {
@@ -784,6 +794,13 @@ void AliTPCclustererMI::Digits2Clusters()
   }  
  
   Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
+
+  if (fUseHLTClusters == 2 && nclusters == 0) {
+    AliInfo("No clusters from TPC Raw data, now trying to read HLT clusters.");
+
+    fZWidth = fParam->GetZWidth();
+    ReadHLTClusters();
+  }
 }
 
 void AliTPCclustererMI::ProcessSectorData(){
@@ -939,6 +956,38 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
   }
   fRowDig = NULL;
 
+  //-----------------------------------------------------------------
+  // Use HLT clusters
+  //-----------------------------------------------------------------
+  if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
+    AliInfo("Using HLT clusters for TPC off-line reconstruction");
+    fZWidth = fParam->GetZWidth();
+    Int_t iResult = ReadHLTClusters();
+
+    // HLT clusters present
+    if (iResult >= 0 && fNclusters > 0)
+      return;
+
+    // HLT clusters not present
+    if (iResult < 0 || fNclusters == 0) {
+      if (fUseHLTClusters == 3) { 
+       AliError("No HLT clusters present, but requested.");
+       return;
+      }
+      else {
+       AliInfo("Now trying to read TPC RAW");
+      }
+    }
+    // Some other problem during cluster reading
+    else {
+      AliWarning("Some problem while unpacking of HLT clusters.");
+      return;
+    }
+  } // if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
+   
+  //-----------------------------------------------------------------
+  // Run TPC off-line clusterer
+  //-----------------------------------------------------------------
   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
   AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
   //
@@ -1106,6 +1155,13 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
   if(fBClonesArray) {
     //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
   }
+
+  if (fUseHLTClusters == 2 && fNclusters == 0) {
+    AliInfo("No clusters from TPC Raw data, now trying to read HLT clusters.");
+
+    fZWidth = fParam->GetZWidth();
+    ReadHLTClusters();
+  }
 }
 
 
@@ -1349,7 +1405,7 @@ void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
 }
 
 Bool_t AliTPCclustererMI::AcceptCluster(AliTPCclusterMI *cl){
-  //
+  // -- Depricated --
   // Currently hack to filter digital noise (15.06.2008)
   // To be parameterized in the AliTPCrecoParam
   // More inteligent way  to be used in future
@@ -1528,7 +1584,111 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t
   return median;
 }
 
+Int_t AliTPCclustererMI::ReadHLTClusters()
+{
+  //
+  // read HLT clusters instead of off line custers, 
+  // used in Digits2Clusters
+  //
 
+  if (!fHLTClusterAccess) {
+  TClass* pCl=NULL;
+  ROOT::NewFunc_t pNewFunc=NULL;
+  do {
+    pCl=TClass::GetClass("AliHLTTPCClusterAccessHLTOUT");
+  } while (!pCl && gSystem->Load("libAliHLTTPC.so")==0);
+  if (!pCl || (pNewFunc=pCl->GetNew())==NULL) {
+    AliError("can not load class description of AliHLTTPCClusterAccessHLTOUT, aborting ...");
+    return -1;
+  }
+  
+  void* p=(*pNewFunc)(NULL);
+  if (!p) {
+    AliError("unable to create instance of AliHLTTPCClusterAccessHLTOUT");
+    return -2;
+  }
+  fHLTClusterAccess=reinterpret_cast<TObject*>(p);
+  }
+
+  TObject* pClusterAccess=fHLTClusterAccess;
+
+  const Int_t kNIS = fParam->GetNInnerSector();
+  const Int_t kNOS = fParam->GetNOuterSector();
+  const Int_t kNS = kNIS + kNOS;
+  fNclusters  = 0;
+  
+  // make sure that all clusters from the previous event are cleared
+  pClusterAccess->Clear("event");
+  for(fSector = 0; fSector < kNS; fSector++) {
+
+    Int_t iResult = 1;
+    TString param("sector="); param+=fSector;
+    // prepare for next sector
+    pClusterAccess->Clear("sector");
+    pClusterAccess->Execute("read", param, &iResult);
+    if (iResult < 0) {
+      return iResult;
+      AliError("HLT Clusters can not be found");
+    }
+
+    TObject* pObj=pClusterAccess->FindObject("clusterarray");
+    if (pObj==NULL) {
+      AliError("HLT clusters requested, but not cluster array not present");
+      return -4;
+    }
 
+    TObjArray* clusterArray=dynamic_cast<TClonesArray*>(pObj);
+    if (!clusterArray) {
+      AliError("HLT cluster array is not of class type TClonesArray");
+      return -5;
+    }
+
+    AliDebug(4,Form("Reading %d clusters from HLT for sector %d", clusterArray->GetEntriesFast(), fSector));
+
+    Int_t nClusterSector=0;
+    Int_t nRows=fParam->GetNRow(fSector);
+
+    for (fRow = 0; fRow < nRows; fRow++) {
+      fRowCl->SetID(fParam->GetIndex(fSector, fRow));
+      if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
+      fNcluster=0; // reset clusters per row
+      
+      fRx = fParam->GetPadRowRadii(fSector, fRow);
+      fPadLength = fParam->GetPadPitchLength(fSector, fRow);
+      fPadWidth  = fParam->GetPadPitchWidth();
+      fMaxPad = fParam->GetNPads(fSector,fRow);
+      fMaxBin = fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
+      
+      fBins = fAllBins[fRow];
+      fSigBins = fAllSigBins[fRow];
+      fNSigBins = fAllNSigBins[fRow];
+      
+      for (Int_t i=0; i<clusterArray->GetEntriesFast(); i++) {
+       if (!clusterArray->At(i)) 
+         continue;
+       
+       AliTPCclusterMI* cluster=dynamic_cast<AliTPCclusterMI*>(clusterArray->At(i));
+       if (!cluster) continue;
+       if (cluster->GetRow()!=fRow) continue;
+       nClusterSector++;
+       AddCluster(*cluster, NULL, 0);
+      }
+      
+      FillRow();
+      fRowCl->GetArray()->Clear("c");
+      
+    } // for (fRow = 0; fRow < nRows; fRow++) {
+    if (nClusterSector!=clusterArray->GetEntriesFast()) {
+      AliError(Form("Failed to read %d out of %d HLT clusters", 
+                   clusterArray->GetEntriesFast()-nClusterSector, 
+                   clusterArray->GetEntriesFast()));
+    }
+    fNclusters+=nClusterSector;
+  } // for(fSector = 0; fSector < kNS; fSector++) {
 
+  pClusterAccess->Clear("event");
 
+  Info("Digits2Clusters", "Number of converted HLT clusters : %d", fNclusters);
+  
+  return 0;
+}