]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCclustererMI.cxx
cout removed
[u/mrichter/AliRoot.git] / TPC / AliTPCclustererMI.cxx
index 3af4a8c422680efa3f4367808a144168482ea1a6..8b481f7b859fcd6a1b14c4fe3753599f10582d6d 100644 (file)
@@ -82,6 +82,8 @@
 #include "AliTPCTransform.h"
 #include "AliTPCclustererMI.h"
 
+using std::cerr;
+using std::endl;
 ClassImp(AliTPCclustererMI)
 
 
@@ -118,10 +120,11 @@ AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoPar
   fRecoParam(0),
   fBDumpSignal(kFALSE),
   fBClonesArray(kFALSE),
-  fUseHLTClusters(1),
+  fUseHLTClusters(4),
   fAllBins(NULL),
   fAllSigBins(NULL),
-  fAllNSigBins(NULL)
+  fAllNSigBins(NULL),
+  fHLTClusterAccess(NULL)
 {
   //
   // COSNTRUCTOR
@@ -165,59 +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),
-  fUseHLTClusters(1),
-  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(){
   //
@@ -247,6 +198,7 @@ AliTPCclustererMI::~AliTPCclustererMI(){
   delete [] fAllBins;
   delete [] fAllSigBins;
   delete [] fAllNSigBins;
+  if (fHLTClusterAccess) delete fHLTClusterAccess;
 }
 
 void AliTPCclustererMI::SetInput(TTree * tree)
@@ -319,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
 
@@ -349,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
@@ -659,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;
@@ -721,34 +686,36 @@ void AliTPCclustererMI::Digits2Clusters()
     fParam->Dump();
     fRecoParam->Dump();
   }
+  fRowDig = NULL;
 
   //-----------------------------------------------------------------
   // Use HLT clusters
   //-----------------------------------------------------------------
-  fUseHLTClusters = fRecoParam->GetUseHLTClusters();
-
-  printf(" HLT TPC Reco foo : %d \n",fUseHLTClusters );
-  if (fUseHLTClusters == 3 && fUseHLTClusters == 4) {
+  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)
-      return;
+    if (iResult >= 0 && fNclusters > 0)
+      return; 
+
     // HLT clusters not present
-    else if(iResult == -1) {
-      if (fUseHLTClusters == 3) {
-       AliError("No HLT clusters present, but requiered.");
+    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) {
+  } // if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
 
   //-----------------------------------------------------------------
   // Run TPC off-line clusterer
@@ -784,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 {
@@ -827,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(){
@@ -985,29 +959,31 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
   //-----------------------------------------------------------------
   // Use HLT clusters
   //-----------------------------------------------------------------
-  fUseHLTClusters = fRecoParam->GetUseHLTClusters();
-
-  if (fUseHLTClusters == 3 && fUseHLTClusters == 4) {
+  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)
+    if (iResult >= 0 && fNclusters > 0)
       return;
+
     // HLT clusters not present
-    else if(iResult == -1) {
-      if (fUseHLTClusters == 3) {
-       AliError("No HLT clusters present, but requiered.");
+    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) {
+  } // if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
    
   //-----------------------------------------------------------------
   // Run TPC off-line clusterer
@@ -1179,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();
+  }
 }
 
 
@@ -1422,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
@@ -1608,7 +1591,7 @@ Int_t AliTPCclustererMI::ReadHLTClusters()
   // used in Digits2Clusters
   //
 
-  TObject* pClusterAccess=NULL;
+  if (!fHLTClusterAccess) {
   TClass* pCl=NULL;
   ROOT::NewFunc_t pNewFunc=NULL;
   do {
@@ -1624,28 +1607,37 @@ Int_t AliTPCclustererMI::ReadHLTClusters()
     AliError("unable to create instance of AliHLTTPCClusterAccessHLTOUT");
     return -2;
   }
-  pClusterAccess=reinterpret_cast<TObject*>(p);
-  if (!pClusterAccess) {
-    AliError("instance not of type TObject");
-    return -3 ;
+  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;
-    pClusterAccess->Clear();
-    pClusterAccess->Execute("read", param);
-    if (pClusterAccess->FindObject("clusterarray")==NULL) {
+    // 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;
     }
 
-    TClonesArray* clusterArray=dynamic_cast<TClonesArray*>(pClusterAccess->FindObject("clusterarray"));
+    TObjArray* clusterArray=dynamic_cast<TClonesArray*>(pObj);
     if (!clusterArray) {
       AliError("HLT cluster array is not of class type TClonesArray");
       return -5;
@@ -1694,7 +1686,7 @@ Int_t AliTPCclustererMI::ReadHLTClusters()
     fNclusters+=nClusterSector;
   } // for(fSector = 0; fSector < kNS; fSector++) {
 
-  delete pClusterAccess;
+  pClusterAccess->Clear("event");
 
   Info("Digits2Clusters", "Number of converted HLT clusters : %d", fNclusters);