]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - HLT/src/AliL3ClustFinderNew.cxx
Solved bug to not merge clusters on the same pad.
[u/mrichter/AliRoot.git] / HLT / src / AliL3ClustFinderNew.cxx
index e05fb823a8f4723bad9f82f9338ca57fa15fb77b..eeb96265cb69fb59ebf30f3b9a16f12b03d361c6 100644 (file)
@@ -1,36 +1,42 @@
-//Author:        Anders Strand Vestbo
-//Last Modified: 17.08.2001
+//$Id$
+
+// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
+//*-- Copyright &copy ASV 
 
 #include <iostream.h>
+#include <math.h>
 #include "AliL3Logging.h"
 #include "AliL3ClustFinderNew.h"
 #include "AliL3DigitData.h"
 #include "AliL3Transform.h"
 #include "AliL3SpacePointData.h"
+#include "AliL3MemHandler.h"
+
+/** \class AliL3ClustFinderNew
+//<pre>
+//_____________________________________________________________
+// AliL3ClustFinderNew
+//
+// The current cluster finder for HLT
+// Based on STAR L3
+//</pre> */
 
 ClassImp(AliL3ClustFinderNew)
 
 AliL3ClustFinderNew::AliL3ClustFinderNew()
 {
   fMatch = 4;
-  fThreshold =10;
+  fThreshold = 10;
+  fXYErr = 0.2;
+  fZErr = 0.3;
   fDeconvPad = kTRUE;
-  fDeconvTime = fTRUE;
-}
-
-AliL3ClustFinderNew::AliL3ClustFinderNew(AliL3Transform *transform)
-{
-  fTransform = transform;
-  fMatch = 4;
-  fThreshold =10;
   fDeconvTime = kTRUE;
-  fDeconvPad = kTRUE;
+  fstdout = kFALSE;
+  fcalcerr = kTRUE;
 }
 
 AliL3ClustFinderNew::~AliL3ClustFinderNew()
 {
-
-
 }
 
 void AliL3ClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_t lastrow,Int_t nmaxpoints)
@@ -41,8 +47,6 @@ void AliL3ClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t firstrow, Int_
   fCurrentPatch = patch;
   fFirstRow = firstrow;
   fLastRow = lastrow;
-  fDeconvTime = kTRUE;
-  fDeconvPad = kTRUE;
 }
 
 void AliL3ClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
@@ -55,16 +59,13 @@ void AliL3ClustFinderNew::InitSlice(Int_t slice,Int_t patch,Int_t nmaxpoints)
 
 void AliL3ClustFinderNew::SetOutputArray(AliL3SpacePointData *pt)
 {
-  
   fSpacePointData = pt;
 }
 
-
 void AliL3ClustFinderNew::Read(UInt_t ndigits,AliL3DigitRowData *ptr)
 {
   fNDigitRowData = ndigits;
   fDigitRowData = ptr;
-  
 }
 
 void AliL3ClustFinderNew::ProcessDigits()
@@ -83,10 +84,11 @@ void AliL3ClustFinderNew::ProcessDigits()
       tmp += size;
       tempPt = (AliL3DigitRowData*)tmp;
     }
+  LOG(AliL3Log::kInformational,"AliL3ClustFinderNew::WriteClusters","Space points")
+    <<"Cluster finder found "<<fNClusters<<" clusters in slice "<<fCurrentSlice
+    <<" patch "<<fCurrentPatch<<ENDLOG;
 }
 
-
-
 void AliL3ClustFinderNew::ProcessRow(AliL3DigitRowData *tempPt)
 {
 
@@ -96,7 +98,7 @@ void AliL3ClustFinderNew::ProcessRow(AliL3DigitRowData *tempPt)
   ClusterData *pad2[2500]; //2 lists for internal memory=2pads
   ClusterData clusterlist[5000]; //Clusterlist
 
-  ClusterData **currentPt; //List of pointers to the current pad
+  ClusterData **currentPt;  //List of pointers to the current pad
   ClusterData **previousPt; //List of pointers to the previous pad
   currentPt = pad2;
   previousPt = pad1;
@@ -105,7 +107,6 @@ void AliL3ClustFinderNew::ProcessRow(AliL3DigitRowData *tempPt)
   //Loop over sequences of this row:
   for(UInt_t bin=0; bin<tempPt->fNDigit; bin++)
     {
-      
       AliL3DigitData *data = tempPt->fDigitData;
       if(data[bin].fPad != last_pad)
        {
@@ -116,7 +117,6 @@ void AliL3ClustFinderNew::ProcessRow(AliL3DigitRowData *tempPt)
            {
              currentPt = pad1;
              previousPt = pad2;
-             
            }
          else 
            {
@@ -134,11 +134,10 @@ void AliL3ClustFinderNew::ProcessRow(AliL3DigitRowData *tempPt)
        }
 
       Bool_t new_cluster = kTRUE;
-
-      UInt_t seq_charge=0,seq_average=0;
-      
+      UInt_t seq_charge=0,seq_average=0,seq_error=0;
       UInt_t last_charge=0,last_was_falling=0;
       Int_t new_bin=-1;
+
       if(fDeconvTime)
        {
        redo: //This is a goto.
@@ -152,9 +151,15 @@ void AliL3ClustFinderNew::ProcessRow(AliL3DigitRowData *tempPt)
          last_was_falling = 0;
        }
       
-
       while(1) //Loop over current sequence
        {
+         if(data[bin].fTime >= AliL3Transform::GetNTimeBins())
+           {
+             LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Digits")
+               <<"Timebin out of range "<<(Int_t)data[bin].fTime<<ENDLOG;
+             break;
+           }
+
          //Get the current ADC-value
          UInt_t charge = data[bin].fCharge;
          
@@ -176,7 +181,8 @@ void AliL3ClustFinderNew::ProcessRow(AliL3DigitRowData *tempPt)
          //Sum the total charge of this sequence
          seq_charge += charge;
          seq_average += data[bin].fTime*charge;
-         
+         seq_error += data[bin].fTime*data[bin].fTime*charge;
+
          //Check where to stop:
          if(data[bin+1].fPad != data[bin].fPad) //new pad
            break; 
@@ -192,32 +198,36 @@ void AliL3ClustFinderNew::ProcessRow(AliL3DigitRowData *tempPt)
        seq_mean = seq_average/seq_charge;
       else
        {
-         printf("\nCF: Should not happen\n");
+         LOG(AliL3Log::kFatal,"AliL3ClustFinderNew::ProcessRow","Data")
+           <<"Error in data given to the cluster finder"<<ENDLOG;
          seq_mean = 1;
          seq_charge = 1;
        }
       
       //Calculate mean in pad direction:
       Int_t pad_mean = seq_charge*data[bin].fPad;
-
+      Int_t pad_error = data[bin].fPad*pad_mean;
+      
       //Compare with results on previous pad:
       for(UInt_t p=0; p<n_previous; p++)
        {
+         //dont merge sequences on the same pad twice
+         if(previousPt[p]->fLastMergedPad==data[bin].fPad) continue;
+
          Int_t difference = seq_mean - previousPt[p]->fMean;
-         if(difference < -fMatch)
-           break;
+         if(difference < -fMatch) break;
 
-         if(difference <= fMatch)//There is a match here!!
+         if(difference <= fMatch) //There is a match here!!
            {
-             
              ClusterData *local = previousPt[p];
+
              if(fDeconvPad)
                {
                  if(seq_charge > local->fLastCharge)
                    {
-                     if(local->fChargeFalling)//The previous pad was falling
+                     if(local->fChargeFalling) //The previous pad was falling
                        {                       
-                         break;//create a new cluster
+                         break; //create a new cluster
                        }                   
                    }
                  else
@@ -229,20 +239,21 @@ void AliL3ClustFinderNew::ProcessRow(AliL3DigitRowData *tempPt)
              new_cluster = kFALSE;
              
              //Update cluster on current pad with the matching one:
-             //currentPt[n_current] = previousPt[p];
-             
              local->fTotalCharge += seq_charge;
              local->fPad += pad_mean;
+             local->fPad2 += pad_error;
              local->fTime += seq_average;
+             local->fTime2 += seq_error;
              local->fMean = seq_mean;
              local->fFlags++; //means we have more than one pad 
-             
+             local->fLastMergedPad = data[bin].fPad;
+
              currentPt[n_current] = local;
              n_current++;
              
              break;
-           }//Checking for match at previous pad
-       }//Loop over results on previous pad.
+           } //Checking for match at previous pad
+       } //Loop over results on previous pad.
       
       if(new_cluster)
        {
@@ -254,9 +265,12 @@ void AliL3ClustFinderNew::ProcessRow(AliL3DigitRowData *tempPt)
          ClusterData *tmp = &clusterlist[n_total];
          tmp->fTotalCharge = seq_charge;
          tmp->fPad = pad_mean;
+         tmp->fPad2 = pad_error;
          tmp->fTime = seq_average;
+         tmp->fTime2 = seq_error;
          tmp->fMean = seq_mean;
          tmp->fFlags = 0;  //flags for 1 pad clusters
+         tmp->fLastMergedPad = data[bin].fPad;
          if(fDeconvPad)
            {
              tmp->fChargeFalling = 0;
@@ -268,14 +282,12 @@ void AliL3ClustFinderNew::ProcessRow(AliL3DigitRowData *tempPt)
          n_total++;
          n_current++;
        }
+
       if(fDeconvTime)
        if(new_bin >= 0) goto redo;
-      
     }//Loop over digits on this padrow
   
   WriteClusters(n_total,clusterlist);
-
-
 }
 
 void AliL3ClustFinderNew::WriteClusters(Int_t n_clusters,ClusterData *list)
@@ -287,13 +299,25 @@ void AliL3ClustFinderNew::WriteClusters(Int_t n_clusters,ClusterData *list)
     {
       if(!list[j].fFlags) continue; //discard 1 pad clusters
       if(list[j].fTotalCharge < fThreshold) continue; //noise cluster
-      Float_t xyz[3];
-      
-      Float_t fpad=(Float_t)list[j].fPad/(Float_t)list[j].fTotalCharge;
-      Float_t ftime=(Float_t)list[j].fTime/(Float_t)list[j].fTotalCharge;
-      //printf("padrow %d number of pads %d totalcharge %d\n",fCurrentRow,list[j].fFlags,list[j].fTotalCharge);
-      fTransform->Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
-      fTransform->Raw2Local(xyz,thissector,thisrow,fpad,ftime);
+
+      Float_t xyz[3];      
+      Float_t fpad =(Float_t)list[j].fPad /(Float_t)list[j].fTotalCharge;
+      Float_t fpad2=fXYErr;
+      Float_t ftime =(Float_t)list[j].fTime /(Float_t)list[j].fTotalCharge;
+      Float_t ftime2=fZErr;
+
+      if(fcalcerr) {
+       fpad2=(Float_t)list[j].fPad2/(Float_t)list[j].fTotalCharge - fpad*fpad;
+       fpad2 = sqrt(fpad2);
+       ftime2=(Float_t)list[j].fTime2/(Float_t)list[j].fTotalCharge - ftime*ftime;
+       ftime2 = sqrt(ftime2); 
+      }
+       
+      if(fstdout==kTRUE)
+       cout<<"WriteCluster: padrow "<<fCurrentRow<<" pad "<<fpad << " +- "<<fpad2<<" time "<<ftime<<" +- "<<ftime2<<" charge "<<list[j].fTotalCharge<<endl;
+
+      AliL3Transform::Slice2Sector(fCurrentSlice,fCurrentRow,thissector,thisrow);
+      AliL3Transform::Raw2Local(xyz,thissector,thisrow,fpad,ftime);
       if(xyz[0]==0) LOG(AliL3Log::kError,"AliL3ClustFinder","Cluster Finder")
        <<AliL3Log::kDec<<"Zero cluster"<<ENDLOG;
       if(fNClusters >= fMaxNClusters)
@@ -301,19 +325,60 @@ void AliL3ClustFinderNew::WriteClusters(Int_t n_clusters,ClusterData *list)
          LOG(AliL3Log::kError,"AliL3ClustFinder::WriteClusters","Cluster Finder")
            <<AliL3Log::kDec<<"Too many clusters"<<ENDLOG;
          return;
-         }  
+       }  
+      fSpacePointData[counter].fCharge = list[j].fTotalCharge;
       fSpacePointData[counter].fX = xyz[0];
       fSpacePointData[counter].fY = xyz[1];
       fSpacePointData[counter].fZ = xyz[2];
       fSpacePointData[counter].fPadRow = fCurrentRow;
-      fSpacePointData[counter].fXYErr = fXYErr;
-      fSpacePointData[counter].fZErr = fZErr;
+      fSpacePointData[counter].fXYErr = fpad2;
+      fSpacePointData[counter].fZErr  = ftime2;
       fSpacePointData[counter].fID = counter
-       +((fCurrentSlice&0x7f)<<25)+((fCurrentPatch&0x7)<<22);//uli
-      
+       +((fCurrentSlice&0x7f)<<25)+((fCurrentPatch&0x7)<<22);//Uli
+#ifdef do_mc
+      Int_t trackID[3];
+      GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
+      //cout<<"padrow "<<fCurrentRow<<" pad "<<(Int_t)rint(fpad)<<" time "<<(Int_t)rint(ftime)<<" Trackid "<<trackID[0]<<endl;
+      fSpacePointData[counter].fTrackID[0] = trackID[0];
+      fSpacePointData[counter].fTrackID[1] = trackID[1];
+      fSpacePointData[counter].fTrackID[2] = trackID[2];
+#endif
       fNClusters++;
       counter++;
-
     }
+}
 
+#ifdef do_mc
+void AliL3ClustFinderNew::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
+{
+  AliL3DigitRowData *rowPt = (AliL3DigitRowData*)fDigitRowData;
+  
+  trackID[0]=trackID[1]=trackID[2]=-2;
+  //cout<<"Looking for pad "<<pad<<" time "<<time<<endl;
+  for(Int_t i=fFirstRow; i<=fLastRow; i++)
+    {
+      if(rowPt->fRow < (UInt_t)fCurrentRow)
+       {
+         AliL3MemHandler::UpdateRowPointer(rowPt);
+         continue;
+       }
+      AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
+      for(UInt_t j=0; j<rowPt->fNDigit; j++)
+       {
+         Int_t cpad = digPt[j].fPad;
+         Int_t ctime = digPt[j].fTime;
+         if(cpad != pad) continue;
+         if(ctime != time) continue;
+         //if(cpad != pad && ctime != ctime) continue;
+         //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
+         trackID[0] = digPt[j].fTrackID[0];
+         trackID[1] = digPt[j].fTrackID[1];
+         trackID[2] = digPt[j].fTrackID[2];
+         break;
+         //cout<<"Reading trackID "<<trackID[0]<<endl;
+       }
+      break;
+    }
+  
 }
+#endif