]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MFT/AliMFTClusterFinder.cxx
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / MFT / AliMFTClusterFinder.cxx
index 165ac720320b31f707d6311a4fc2b274dc057655..f09aa007bbcf426e775c4f93b610b901f4500d9c 100644 (file)
 #include "AliMFTSegmentation.h"
 #include "TTree.h"
 #include "TMath.h"
+#include "AliMFTConstants.h"
 #include "AliMFTClusterFinder.h"
+#include "TRandom.h"
+
+const Double_t AliMFTClusterFinder::fCutForAvailableDigits = AliMFTConstants::fCutForAvailableDigits;
+const Double_t AliMFTClusterFinder::fCutForAttachingDigits = AliMFTConstants::fCutForAttachingDigits;
+const Double_t AliMFTClusterFinder::fMisalignmentMagnitude = AliMFTConstants::fMisalignmentMagnitude;
 
 ClassImp(AliMFTClusterFinder)
 
@@ -38,28 +44,56 @@ ClassImp(AliMFTClusterFinder)
 AliMFTClusterFinder::AliMFTClusterFinder() : 
   TObject(),
   fDigitsInCluster(0),
-  fCurrentDig(0),
+  fCurrentDigit(0),
+  fCurrentCluster(0),
   fSegmentation(0),
-  fNPlanes(0)
+  fNPlanes(0),
+  fApplyMisalignment(kFALSE),
+  fStopWatch(0)
 {
 
   // Default constructor
   
   for (Int_t iPlane=0; iPlane<fNMaxPlanes; iPlane++) fClustersPerPlane[iPlane] = NULL;
   fDigitsInCluster = new TClonesArray("AliMFTDigit", fNMaxDigitsPerCluster);
+  fDigitsInCluster -> SetOwner(kTRUE);
+  fStopWatch = new TStopwatch();
+  fStopWatch -> Reset();
 
 }
 
 //====================================================================================================================================================
 
 AliMFTClusterFinder::~AliMFTClusterFinder() {
-
+  
   AliDebug(1, "Deleting AliMFTClusterFinder...");
+  
+  for (Int_t iPlane=0; iPlane<fNMaxPlanes; iPlane++) {
+    if (fClustersPerPlane[iPlane]) fClustersPerPlane[iPlane]->Delete(); delete fClustersPerPlane[iPlane]; fClustersPerPlane[iPlane] = 0x0;
+  }
+  if (fDigitsInCluster) fDigitsInCluster->Delete(); delete fDigitsInCluster; fDigitsInCluster = NULL;
+  delete fSegmentation; fSegmentation=NULL;
 
-  for (Int_t iPlane=0; iPlane<fNMaxPlanes; iPlane++) delete fClustersPerPlane[iPlane];
-
+  delete fStopWatch;
+  
   AliDebug(1, "... done!");
+  
+}
+
+//====================================================================================================================================================
 
+void AliMFTClusterFinder::Clear(const Option_t* /*opt*/) {
+  
+  AliDebug(1, "Clearing AliMFTClusterFinder...");
+  
+  for (Int_t iPlane=0; iPlane<fNMaxPlanes; iPlane++) {
+    if(fClustersPerPlane[iPlane]) fClustersPerPlane[iPlane]->Clear("C"); 
+  }
+  if(fDigitsInCluster) fDigitsInCluster->Clear("C");   
+  fSegmentation->Clear("C"); 
+  
+  AliDebug(1, "... done!");
+  
 }
 
 //====================================================================================================================================================
@@ -80,7 +114,7 @@ void AliMFTClusterFinder::StartEvent() {
   AliDebug(1, "Starting Event...");
   
   for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
-    fClustersPerPlane[iPlane]->Clear();
+    fClustersPerPlane[iPlane]->Delete();
   }
 
   AliDebug(1, "... done!");
@@ -96,120 +130,126 @@ void AliMFTClusterFinder::DigitsToClusters(const TObjArray *pDigitList) {
   AliInfo("Starting Clusterization for MFT");
   AliDebug(1, Form("nPlanes = %d", fNPlanes));
 
+  if (!fStopWatch) fStopWatch = new TStopwatch();
+  fStopWatch -> Reset();
+
   StartEvent(); 
+  Bool_t isDigAvailableForNewCluster = kTRUE;
   
+  TClonesArray *myDigitList = 0;
+
   for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
 
     AliDebug(1, Form("Plane %02d", iPlane));
 
-    TClonesArray *myDigitList = (TClonesArray*) pDigitList->At(iPlane);
-    //    myDigitList->Sort();
-
-    AliDebug(1, Form("myDigitList->GetEntries() = %d", myDigitList->GetEntries()));
-
-    while (myDigitList->GetEntries()) {
-
-      fDigitsInCluster->Clear();
-      
-      Bool_t clusterUpdated=kTRUE;
-
-      //---------------------------------------------------------------------------------------------------------
-
-      while (clusterUpdated) {    // repeat the loop on the digits until no new compatible digit is found
-
-       clusterUpdated = kFALSE;
-
-       for (Int_t iDig=0; iDig<myDigitList->GetEntries(); iDig++) {
-         fCurrentDig = (AliMFTDigit*) myDigitList->At(iDig);
-         if (fDigitsInCluster->GetEntries()<fNMaxDigitsPerCluster) {
-           if (fDigitsInCluster->GetEntries()==0) {
-             new ((*fDigitsInCluster)[fDigitsInCluster->GetEntries()]) AliMFTDigit(*fCurrentDig);
-             myDigitList->Remove(fCurrentDig);
-             clusterUpdated = kTRUE;
-           }
-           else if (IsCurrentDigitCompatible()) {
-             new ((*fDigitsInCluster)[fDigitsInCluster->GetEntries()]) AliMFTDigit(*fCurrentDig);
-             myDigitList->Remove(fCurrentDig);
-             clusterUpdated = kTRUE;
-           }
-         }
-       }
-       
-       if (clusterUpdated) myDigitList->Compress();
-
-      }
-
-      //---------------------------------------------------------------------------------------------------------
-
-      AliDebug(1, "Building new cluster");
+    Int_t nDetElem = fSegmentation->GetPlane(iPlane)->GetNActiveElements();
+    TClonesArray *clustersPerDetElem[fNMaxDetElemPerPlane] = {0};
+    for (Int_t iDetElem=0; iDetElem<nDetElem; iDetElem++) clustersPerDetElem[iDetElem] = new TClonesArray("AliMFTCluster");
 
-      BuildNewCluster(iPlane);  // here the new cluster is built
+    myDigitList = (TClonesArray*) pDigitList->At(iPlane);
 
-    }
+    AliDebug(1, Form("myDigitList->GetEntries() = %d", myDigitList->GetEntries()));
 
-  }
+    Int_t cycleOverDigits = 0;
+    Double_t myCutForAvailableDigits = 0;
+    
+    Int_t currentDetElem = -1;
+    Int_t currentDetElemLocal = -1;
+    Bool_t areThereSkippedDigits = kFALSE;
 
-}
+    fStopWatch -> Start();
 
-//====================================================================================================================================================
-
-Bool_t AliMFTClusterFinder::IsCurrentDigitCompatible() {
+    while (myDigitList->GetEntries()) {
 
-  // where it is decided if the current digit (fCurrentDig) is compatible with the current digits array (fDigitsInCluster)
+      for (Int_t iDig=0; iDig<myDigitList->GetEntries(); iDig++) {
 
-  for (Int_t iDigit=0; iDigit<fDigitsInCluster->GetEntries(); iDigit++) {
-    AliMFTDigit *tmpDig = (AliMFTDigit*) fDigitsInCluster->At(iDigit);
-    Int_t distX = TMath::Abs(tmpDig->GetPixelX() - fCurrentDig->GetPixelX());
-    Int_t distY = TMath::Abs(tmpDig->GetPixelY() - fCurrentDig->GetPixelY());
-    if (distX<=1 &&  distY<=1) return kTRUE;
-  }
+       fCurrentDigit = (AliMFTDigit*) myDigitList->At(iDig);
 
-  return kFALSE;
-    
-}
+       if (!iDig) {
+         if (fCurrentDigit->GetDetElemID() != currentDetElem) {      
+           // first iteration over the digits of a specific detection element
+           currentDetElem = fCurrentDigit->GetDetElemID();
+           currentDetElemLocal = fSegmentation->GetDetElemLocalID(currentDetElem);
+           cycleOverDigits = 0;
+           myCutForAvailableDigits = fCutForAvailableDigits;
+         }
+         else if (fCurrentDigit->GetDetElemID()==currentDetElem && areThereSkippedDigits) {
+           // second (or further) iteration over the digits of a specific detection element
+           cycleOverDigits++;
+           myCutForAvailableDigits -= 0.5;
+         }
+         areThereSkippedDigits = kFALSE;
+       }
+       else {
+         areThereSkippedDigits = kTRUE;
+         if (fCurrentDigit->GetDetElemID() != currentDetElem) break;
+       }
 
-//====================================================================================================================================================
+       isDigAvailableForNewCluster = kTRUE;
+
+       for (Int_t iCluster=0; iCluster<clustersPerDetElem[currentDetElemLocal]->GetEntries(); iCluster++) {
+         fCurrentCluster = (AliMFTCluster*) clustersPerDetElem[currentDetElemLocal]->At(iCluster);
+         if (fCurrentCluster->GetDistanceFromPixel(fCurrentDigit) < fCutForAttachingDigits) { 
+           fCurrentCluster->AddPixel(fCurrentDigit); 
+           myDigitList->Remove(fCurrentDigit);
+           myDigitList->Compress();
+           iDig--;
+           isDigAvailableForNewCluster = kFALSE;
+           break; 
+         }
+         if (fCurrentCluster->GetDistanceFromPixel(fCurrentDigit) < myCutForAvailableDigits) isDigAvailableForNewCluster=kFALSE;
+       }
 
-void AliMFTClusterFinder::BuildNewCluster(Int_t plane) {
+       if (isDigAvailableForNewCluster) {
+         AliMFTCluster *newCluster = new AliMFTCluster();
+         newCluster->AddPixel(fCurrentDigit);
+         myDigitList->Remove(fCurrentDigit);
+         myDigitList->Compress();
+         iDig--;
+         new ((*clustersPerDetElem[currentDetElemLocal])[clustersPerDetElem[currentDetElemLocal]->GetEntries()]) AliMFTCluster(*newCluster);
+         delete newCluster;
+       }
+       
+      }   // end of cycle over the digits
 
-  // where a new cluster is built, starting from the array of compatible digits (fDigitsInCluster)
+    }   // no more digits to check in current plane!
 
-  AliDebug(1, Form("Starting cluster building from %d digits", fDigitsInCluster->GetEntries()));
+    fStopWatch -> Print("m");
 
-  AliMFTCluster *newCluster = new AliMFTCluster();
+    printf("Plane %d: clusters found in %f seconds\n",iPlane,fStopWatch->CpuTime());
+    fStopWatch->Start();
 
-  Double_t xCenters[fNMaxDigitsPerCluster] = {0};
-  Double_t yCenters[fNMaxDigitsPerCluster] = {0};
-  Double_t nElectrons = 0.;
+    // Now we merge the cluster lists coming from each detection element, to build the cluster list of the plane
 
-  for (Int_t iDigit=0; iDigit<fDigitsInCluster->GetEntries(); iDigit++) {
-    AliMFTDigit *tmpDig = (AliMFTDigit*) fDigitsInCluster->At(iDigit);
-    xCenters[iDigit] = tmpDig->GetPixelCenterX();
-    yCenters[iDigit] = tmpDig->GetPixelCenterY();
-    nElectrons      += tmpDig->GetNElectrons();
-    for (Int_t iTrack=0; iTrack<tmpDig->GetNMCTracks(); iTrack++) newCluster->AddMCLabel(tmpDig->GetMCLabel(iTrack));
-  }
+    Double_t misalignmentPhi = 2.*TMath::Pi()*gRandom->Rndm();
+    Double_t misalignmentX   = fMisalignmentMagnitude*TMath::Cos(misalignmentPhi);
+    Double_t misalignmentY   = fMisalignmentMagnitude*TMath::Sin(misalignmentPhi);
 
-  newCluster -> SetX(TMath::Mean(fDigitsInCluster->GetEntries(), xCenters));
-  newCluster -> SetY(TMath::Mean(fDigitsInCluster->GetEntries(), yCenters));
-  newCluster -> SetZ(((AliMFTDigit*) fDigitsInCluster->At(0))->GetPixelCenterZ());
+    AliMFTCluster *newCluster = NULL;
+    for (Int_t iDetElem=0; iDetElem<nDetElem; iDetElem++) {
+      for (Int_t iCluster=0; iCluster<clustersPerDetElem[iDetElem]->GetEntries(); iCluster++) {
+       newCluster = (AliMFTCluster*) (clustersPerDetElem[iDetElem]->At(iCluster));
+       newCluster -> TerminateCluster();
+       if (fApplyMisalignment) {
+         newCluster -> SetClusterEditable(kTRUE);
+         newCluster -> SetX(newCluster->GetX()+misalignmentX);
+         newCluster -> SetY(newCluster->GetY()+misalignmentY);
+         newCluster -> SetClusterEditable(kFALSE);
+       }
+       new ((*fClustersPerPlane[iPlane])[fClustersPerPlane[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
+      }
+    }
 
-  Double_t minErrX = ((AliMFTDigit*) fDigitsInCluster->At(0))->GetPixelWidthX() / TMath::Sqrt(12.);
-  Double_t minErrY = ((AliMFTDigit*) fDigitsInCluster->At(0))->GetPixelWidthY() / TMath::Sqrt(12.);
-  Double_t minErrZ = ((AliMFTDigit*) fDigitsInCluster->At(0))->GetPixelWidthZ() / TMath::Sqrt(12.);
-  newCluster -> SetErrX( TMath::Max(TMath::RMS(fDigitsInCluster->GetEntries(), xCenters), minErrX) );
-  newCluster -> SetErrY( TMath::Max(TMath::RMS(fDigitsInCluster->GetEntries(), yCenters), minErrY) );
-  newCluster -> SetErrZ( minErrZ );
-    
-  newCluster -> SetNElectrons(nElectrons);
-  newCluster -> SetPlane(plane);
+    printf("%d Clusters in plane %02d merged in %f seconds\n", fClustersPerPlane[iPlane]->GetEntries(), iPlane, fStopWatch->CpuTime());
 
-  newCluster -> SetSize(fDigitsInCluster->GetEntries());
+    for (Int_t iDetElem=0; iDetElem<nDetElem; iDetElem++) {
+      clustersPerDetElem[iDetElem] -> Delete();
+      delete clustersPerDetElem[iDetElem];
+    }
 
-  AliDebug(1, Form("Adding cluster #%02d to tree (%f, %f, %f)", 
-                  fClustersPerPlane[plane]->GetEntries(), newCluster->GetX(), newCluster->GetY(), newCluster->GetZ()));
+    myDigitList -> Delete();
 
-  new ((*fClustersPerPlane[plane])[fClustersPerPlane[plane]->GetEntries()]) AliMFTCluster(*newCluster);
+  }  // end of cycle over the planes
 
 }
 
@@ -226,10 +266,9 @@ void AliMFTClusterFinder::MakeClusterBranch(TTree *treeCluster) {
   if (treeCluster) {
     for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
       AliDebug(1, Form("Setting Branch Plane_%02d for treeCluster",iPlane));
-      TBranch *branch = treeCluster->GetBranch(Form("Plane_%02d",iPlane));
-      if (branch) continue;
+      if (treeCluster->GetBranch(Form("Plane_%02d",iPlane))) continue;
       AliDebug(1, Form("Branch Plane_%02d does not exist, creating!",iPlane));
-      branch = treeCluster->Branch(Form("Plane_%02d",iPlane), &(fClustersPerPlane[iPlane])); 
+      treeCluster->Branch(Form("Plane_%02d",iPlane), &(fClustersPerPlane[iPlane])); 
     }
   }
 
@@ -266,6 +305,8 @@ void AliMFTClusterFinder::CreateClusters() {
   for(Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
     AliDebug(1, Form("plane %02d", iPlane));
     fClustersPerPlane[iPlane] = new TClonesArray("AliMFTCluster");
+    fClustersPerPlane[iPlane] -> SetOwner(kTRUE);
+
   }
 
 }