fix finding of pad neighbours; remove methods to write them in OCDB
[u/mrichter/AliRoot.git] / MUON / AliMUONSimpleClusterServer.cxx
index 542661b..b6bfe99 100644 (file)
 
 #include "AliMUONSimpleClusterServer.h"
 
-#include "AliMUONConstants.h"
+#include "AliCodeTimer.h"
+#include "AliLog.h"
 #include "AliMUONCluster.h"
+#include "AliMUONConstants.h"
 #include "AliMUONGeometryTransformer.h"
 #include "AliMUONPad.h"
+#include "AliMUONTriggerTrackToTrackerClusters.h"
 #include "AliMUONVCluster.h"
 #include "AliMUONVClusterFinder.h"
 #include "AliMUONVClusterStore.h"
 #include "AliMUONVDigitStore.h"
+#include "AliMUONRecoParam.h"
 #include "AliMpArea.h"
 #include "AliMpDEIterator.h"
 #include "AliMpDEManager.h"
 #include "AliMpExMap.h"
+#include "AliMpExMapIterator.h"
+#include "AliMpPad.h"
 #include "AliMpSegmentation.h"
-#include "AliMpVPadIterator.h"
 #include "AliMpVSegmentation.h"
-#include "AliLog.h"
-#include <float.h>
 #include <Riostream.h>
-#include <TClonesArray.h>
+#include <TObjArray.h>
 #include <TString.h>
+#include <float.h>
 
 /// \class AliMUONSimpleClusterServer
 ///
@@ -45,6 +49,8 @@
 /// 
 /// \author Laurent Aphecetche, Subatech
 
+using std::endl;
+using std::cout;
 /// \cond CLASSIMP  
 ClassImp(AliMUONSimpleClusterServer)
 /// \endcond
@@ -54,106 +60,52 @@ namespace
   TString AsString(const AliMpArea& area)
   {
     return Form("(X,Y)=(%7.3f,%7.3f) (DX,DY)=(%7.3f,%7.3f)",
-                area.Position().X(),
-                area.Position().Y(),
-                area.Dimensions().Y(),
-                area.Dimensions().Y());
+                area.GetPositionX(),
+                area.GetPositionY(),
+                area.GetDimensionX(),  /// TBCL was Y !!!
+                area.GetDimensionY());
   }
 }
 
 //_____________________________________________________________________________
-AliMUONSimpleClusterServer::AliMUONSimpleClusterServer(AliMUONVClusterFinder& clusterFinder,
+AliMUONSimpleClusterServer::AliMUONSimpleClusterServer(AliMUONVClusterFinder* clusterFinder,
                                                        const AliMUONGeometryTransformer& transformer)
-: AliMUONVClusterServer(), 
+: AliMUONVClusterServer(),
+  fDigitStore(0x0),
   fClusterFinder(clusterFinder),
-  fTransformer(transformer),
+  fkTransformer(transformer),
   fPads(),
-  fDEAreas(new AliMpExMap(true))
+  fTriggerTrackStore(0x0),
+  fBypass(0x0)
 {
     /// Ctor
+    /// Note that we take ownership of the clusterFinder
     
-    AliMpDEIterator it;
-    
-    it.First();
-    
-    /// Generate the DE areas in global coordinates
+    fPads[0] = new AliMpExMap;
+    fPads[1] = new AliMpExMap;
     
-    while ( !it.IsDone() )
-    {
-      Int_t detElemId = it.CurrentDEId();
-      
-      const AliMpVSegmentation* seg = AliMpSegmentation::Instance()->GetMpSegmentation(detElemId,AliMp::kCath0);
-      
-      Double_t xg,yg,zg;
-      
-      AliMp::StationType stationType = AliMpDEManager::GetStationType(detElemId);
-      
-      Double_t xl(0.0), yl(0.0), zl(0.0);
-      Double_t dx(seg->Dimensions().X());
-      Double_t dy(seg->Dimensions().Y());
-      
-      if ( stationType == AliMp::kStation1 || stationType == AliMp::kStation2 ) 
-      {
-        Double_t xmin(FLT_MAX);
-        Double_t xmax(-FLT_MAX);
-        Double_t ymin(FLT_MAX);
-        Double_t ymax(-FLT_MAX);
-        
-        for ( Int_t icathode = 0; icathode < 2; ++icathode ) 
-        {
-          const AliMpVSegmentation* cathode 
-            = AliMpSegmentation::Instance()->GetMpSegmentation(detElemId,AliMp::GetCathodType(icathode));
-          
-          AliMpVPadIterator* it = cathode->CreateIterator();
-          
-          it->First();
-          
-          while ( !it->IsDone() ) 
-          {
-            AliMpPad pad = it->CurrentItem();
-            AliMpArea a(pad.Position(),pad.Dimensions());
-            xmin = TMath::Min(xmin,a.LeftBorder());
-            xmax = TMath::Max(xmax,a.RightBorder());
-            ymin = TMath::Min(ymin,a.DownBorder());
-            ymax = TMath::Max(ymax,a.UpBorder());
-            it->Next();
-          }
-          
-          delete it;
-        }
-        
-        xl = (xmin+xmax)/2.0;
-        yl = (ymin+ymax)/2.0;
-        dx = (xmax-xmin)/2.0;
-        dy = (ymax-ymin)/2.0;
-        
-        fTransformer.Local2Global(detElemId,xl,yl,zl,xg,yg,zg);
-      }
-      else
-      {
-        fTransformer.Local2Global(detElemId,xl,yl,zl,xg,yg,zg);
-      }
-      
-      fDEAreas->Add(detElemId,new AliMpArea(TVector2(xg,yg),TVector2(dx,dy)));
-     
-      it.Next();
-    }
+    fPadsIterator[0] = fPads[0]->CreateIterator();
+    fPadsIterator[1] = fPads[1]->CreateIterator();
 }
 
 //_____________________________________________________________________________
 AliMUONSimpleClusterServer::~AliMUONSimpleClusterServer()
 {
   /// Dtor
+  delete fClusterFinder;
   delete fPads[0];
   delete fPads[1];
-  delete fDEAreas;
+  delete fPadsIterator[0];
+  delete fPadsIterator[1];
+  delete fBypass;
 }
 
 //_____________________________________________________________________________
 Int_t 
 AliMUONSimpleClusterServer::Clusterize(Int_t chamberId,
                                        AliMUONVClusterStore& clusterStore,
-                                       const AliMpArea& area)
+                                       const AliMpArea& area,
+                                       const AliMUONRecoParam* recoParam)
 {
   /// Area is in absolute coordinate. If not valid, means clusterize all
   /// the chamber.
@@ -161,6 +113,18 @@ AliMUONSimpleClusterServer::Clusterize(Int_t chamberId,
   /// We first find out the list of DE that have a non-zero overlap with area,
   /// and then use the clusterfinder to find clusters in those areas (and DE).
   
+  AliCodeTimerAuto(Form("Chamber %d",chamberId),0);
+  
+  if ( fTriggerTrackStore && chamberId >= 6 ) 
+  {
+    return fBypass->GenerateClusters(chamberId,clusterStore);
+  }
+  
+  if (!recoParam) {
+    AliError("Reconstruction parameters are missing: unable to clusterize");
+    return 0;
+  }
+  
   AliMpDEIterator it;
   
   it.First(chamberId);
@@ -175,10 +139,10 @@ AliMUONSimpleClusterServer::Clusterize(Int_t chamberId,
   {
     Int_t detElemId = it.CurrentDEId();
     
-    TClonesArray* pads[2] = 
+    TObjArray* pads[2] = 
     { 
-      static_cast<TClonesArray*>(fPads[0]->GetValue(detElemId)),
-      static_cast<TClonesArray*>(fPads[1]->GetValue(detElemId)) 
+      static_cast<TObjArray*>(fPads[0]->GetValue(detElemId)),
+      static_cast<TObjArray*>(fPads[1]->GetValue(detElemId)) 
     };
     
     if ( ( pads[0] && pads[0]->GetLast()>=0 ) || 
@@ -194,17 +158,21 @@ AliMUONSimpleClusterServer::Clusterize(Int_t chamberId,
       
       if ( ok ) 
       {      
-        if ( fClusterFinder.NeedSegmentation() )
-        {
-          const AliMpVSegmentation* seg[2] = 
+       const AliMpVSegmentation* seg[2] = 
         { AliMpSegmentation::Instance()->GetMpSegmentation(detElemId,AliMp::kCath0),
           AliMpSegmentation::Instance()->GetMpSegmentation(detElemId,AliMp::kCath1)
         };
-          fClusterFinder.Prepare(detElemId,pads,deArea,seg);
+        
+        fClusterFinder->SetChargeHints(recoParam->LowestPadCharge(),
+                                       recoParam->LowestClusterCharge());
+        
+        if ( fClusterFinder->NeedSegmentation() )
+        {
+          fClusterFinder->Prepare(detElemId,pads,deArea,seg);
         }
         else
         {
-          fClusterFinder.Prepare(detElemId,pads,deArea);
+          fClusterFinder->Prepare(detElemId,pads,deArea);
         }
         
         AliDebug(1,Form("Clusterizing DE %04d with %3d pads (cath0) and %3d pads (cath1)",
@@ -214,11 +182,11 @@ AliMUONSimpleClusterServer::Clusterize(Int_t chamberId,
         
         AliMUONCluster* cluster;
         
-        while ( ( cluster = fClusterFinder.NextCluster() ) ) 
+        while ( ( cluster = fClusterFinder->NextCluster() ) ) 
         {      
           // add new cluster to the store with information to build its ID
           // increment the number of clusters into the store
-          AliMUONVCluster* rawCluster = clusterStore.Add(AliMpDEManager::GetChamberId(detElemId), detElemId, fNCluster++);
+          AliMUONVCluster* rawCluster = clusterStore.Add(chamberId, detElemId, fNCluster++);
           
           ++nofAddedClusters;
           
@@ -229,22 +197,34 @@ AliMUONSimpleClusterServer::Clusterize(Int_t chamberId,
           for (Int_t iPad=0; iPad<nPad; iPad++) 
           {
             AliMUONPad *pad = cluster->Pad(iPad);
-            rawCluster->AddDigitId(pad->GetUniqueID());
+           
+           // skip virtual pads
+           if (!pad->IsReal()) continue;
+            
+           rawCluster->AddDigitId(pad->GetUniqueID());
           }
           
           // fill charge and other cluster informations
           rawCluster->SetCharge(cluster->Charge());
+          rawCluster->SetChi2(cluster->Chi2());
           
           Double_t xg, yg, zg;
-          fTransformer.Local2Global(detElemId, 
+          fkTransformer.Local2Global(detElemId, 
                                     cluster->Position().X(), cluster->Position().Y(), 
                                     0, xg, yg, zg);
           rawCluster->SetXYZ(xg, yg, zg);
+          rawCluster->SetErrXY(recoParam->GetDefaultNonBendingReso(chamberId),recoParam->GetDefaultBendingReso(chamberId));
+          
+          // Set MC label
+          if (fDigitStore && fDigitStore->HasMCInformation()) 
+          {
+            rawCluster->SetMCLabel(FindMCLabel(*cluster, detElemId, seg));
+          }
           
-          AliDebug(1,Form("Adding RawCluster detElemId %4d mult %2d charge %e (xl,yl,zl)=(%e,%e,%e) (xg,yg,zg)=(%e,%e,%e)",
-                          detElemId,nPad,cluster->Charge(),
+          AliDebug(1,Form("Adding RawCluster detElemId %4d mult %2d charge %e (xl,yl,zl)=(%e,%e,%e) (xg,yg,zg)=(%e,%e,%e) label %d",
+                          detElemId,rawCluster->GetNDigits(),rawCluster->GetCharge(),
                           cluster->Position().X(),cluster->Position().Y(),0.0,
-                          xg,yg,zg));
+                          xg,yg,zg,rawCluster->GetMCLabel()));
         }
       }
     }
@@ -252,10 +232,11 @@ AliMUONSimpleClusterServer::Clusterize(Int_t chamberId,
   }
   
   AliDebug(1,Form("chamberId = %2d NofClusters after = %d",chamberId,fNCluster));
-
+  
   return nofAddedClusters;
 }
 
+
 //_____________________________________________________________________________
 void
 AliMUONSimpleClusterServer::Global2Local(Int_t detElemId, const AliMpArea& globalArea,
@@ -265,13 +246,18 @@ AliMUONSimpleClusterServer::Global2Local(Int_t detElemId, const AliMpArea& globa
   
   Double_t xl,yl,zl;
   
-  Double_t zg = AliMUONConstants::DefaultChamberZ(AliMpDEManager::GetChamberId(detElemId));
+  Int_t chamberId = AliMpDEManager::GetChamberId(detElemId);
+  if ( chamberId < 0 ) {
+    AliErrorStream() << "Cannot get chamberId from detElemId=" << detElemId << endl;
+    return;
+  }  
+  Double_t zg = AliMUONConstants::DefaultChamberZ(chamberId);
   
-  fTransformer.Global2Local(detElemId,
-                             globalArea.Position().X(),globalArea.Position().Y(),zg,
+  fkTransformer.Global2Local(detElemId,
+                             globalArea.GetPositionX(),globalArea.GetPositionY(),zg,
                              xl,yl,zl);
   
-  localArea = AliMpArea(TVector2(xl,yl), globalArea.Dimensions());
+  localArea = AliMpArea(xl,yl, globalArea.GetDimensionX(), globalArea.GetDimensionY());
 }
 
 //_____________________________________________________________________________
@@ -286,7 +272,9 @@ AliMUONSimpleClusterServer::Overlap(Int_t detElemId,
   
   Bool_t overlap(kFALSE);
   
-  AliMpArea* globalDEArea = static_cast<AliMpArea*>(fDEAreas->GetValue(detElemId));
+  AliMpArea* globalDEArea = fkTransformer.GetDEArea(detElemId);
+  
+  if (!globalDEArea) return kFALSE;
   
   AliMpArea overlapArea;
   
@@ -313,58 +301,145 @@ AliMUONSimpleClusterServer::Overlap(Int_t detElemId,
 }
 
 //_____________________________________________________________________________
-TClonesArray* 
+TObjArray* 
 AliMUONSimpleClusterServer::PadArray(Int_t detElemId, Int_t cathode) const
 {
   /// Return array for given cathode of given DE
   
-  return static_cast<TClonesArray*>(fPads[cathode]->GetValue(detElemId));
+  return static_cast<TObjArray*>(fPads[cathode]->GetValue(detElemId));
+}
+
+//_____________________________________________________________________________
+Bool_t 
+AliMUONSimpleClusterServer::UseTriggerTrackStore(AliMUONVTriggerTrackStore* trackStore)
+{
+  /// Tells us to use trigger track store, and thus to bypass St45 clusters
+  fTriggerTrackStore = trackStore; // not owner
+  delete fBypass;
+  fBypass = new AliMUONTriggerTrackToTrackerClusters(fkTransformer,fTriggerTrackStore);
+  return kTRUE;
 }
 
 //_____________________________________________________________________________
 void 
-AliMUONSimpleClusterServer::UseDigitStore(const AliMUONVDigitStore& digitStore)
+AliMUONSimpleClusterServer::UseDigits(TIter& next, AliMUONVDigitStore* digitStore)
 {
   /// Convert digitStore into two arrays of AliMUONPads
-
-  delete fPads[0];
-  delete fPads[1];
   
-  fPads[0] = new AliMpExMap(true);
-  fPads[1] = new AliMpExMap(true);
+  fDigitStore = digitStore;
   
-  TIter next(digitStore.CreateIterator());
+  // Clear pads arrays in the maps
+  for ( Int_t i=0; i<2; i++ ) {
+    fPadsIterator[i]->Reset();
+    Int_t key; TObject* obj;
+    while ( ( obj = fPadsIterator[i]->Next(key) ) ) {
+      //cout << "clearing array for detElemId " << key <<  "  ";
+      obj->Clear();
+    }  
+  }   
+
   AliMUONVDigit* d;
-  
   while ( ( d = static_cast<AliMUONVDigit*>(next()) ) )
   {
     if ( ! d->Charge() > 0 ) continue; // skip void digits.
+    if ( ! d->IsTracker() ) continue; // skip trigger digits
     Int_t ix = d->PadX();
     Int_t iy = d->PadY();
     Int_t cathode = d->Cathode();
     Int_t detElemId = d->DetElemId();
     const AliMpVSegmentation* seg = AliMpSegmentation::Instance()->
       GetMpSegmentation(detElemId,AliMp::GetCathodType(cathode));
-    AliMpPad pad = seg->PadByIndices(AliMpIntPair(ix,iy));
+    AliMpPad pad = seg->PadByIndices(ix,iy);
     
-    TClonesArray* padArray = PadArray(detElemId,cathode);
+    TObjArray* padArray = PadArray(detElemId,cathode);
     if (!padArray)
     {
-      padArray = new TClonesArray("AliMUONPad",100);
+      padArray = new TObjArray(100);
+      padArray->SetOwner(kTRUE);
       fPads[cathode]->Add(detElemId,padArray);
     }
     
-    AliMUONPad mpad(detElemId,cathode,
-                    ix,iy,pad.Position().X(),pad.Position().Y(),
-                    pad.Dimensions().X(),pad.Dimensions().Y(),
+    AliMUONPad* mpad = new AliMUONPad(detElemId,cathode,
+                    ix,iy,pad.GetPositionX(),pad.GetPositionY(),
+                    pad.GetDimensionX(),pad.GetDimensionY(),
                     d->Charge());
-    if ( d->IsSaturated() ) mpad.SetSaturated(kTRUE);
-    mpad.SetUniqueID(d->GetUniqueID());
-    new ((*padArray)[padArray->GetLast()+1]) AliMUONPad(mpad);      
+    if ( d->IsSaturated() ) mpad->SetSaturated(kTRUE);
+    mpad->SetUniqueID(d->GetUniqueID());
+    padArray->Add(mpad);      
   }
 }
 
 //_____________________________________________________________________________
+Int_t
+AliMUONSimpleClusterServer::FindMCLabel(const AliMUONCluster& cluster, Int_t detElemId, const AliMpVSegmentation* seg[2]) const
+{
+  /// Find the label of the most contributing MC track (-1 in case of failure)
+  /// The data member fDigitStore must be set
+  
+  // --- get the digit (if any) located at the cluster position on both cathods ---
+  Int_t nTracks[2] = {0, 0};
+  AliMUONVDigit* digit[2] = {0x0, 0x0};
+  for (Int_t iCath = 0; iCath < 2; iCath++) {
+    AliMpPad pad 
+      = seg[AliMp::GetCathodType(iCath)]->PadByPosition(cluster.Position().X(), cluster.Position().Y(),kFALSE);
+    if (pad.IsValid()) {
+      digit[iCath] = fDigitStore->FindObject(detElemId, pad.GetManuId(), pad.GetManuChannel(), iCath);
+      if (digit[iCath]) nTracks[iCath] = digit[iCath]->Ntracks();
+    }
+  }
+  
+  if (nTracks[0] + nTracks[1] == 0) return -1;
+  
+  // --- build the list of contributing tracks and of the associated charge ---
+  Int_t* trackId = new Int_t[nTracks[0] + nTracks[1]];
+  Float_t* trackCharge = new Float_t[nTracks[0] + nTracks[1]];
+  Int_t nTracksTot = 0;
+  
+  // fill with contributing tracks on first cathod
+  for (Int_t iTrack1 = 0; iTrack1 < nTracks[0]; iTrack1++) {
+    trackId[iTrack1] = digit[0]->Track(iTrack1);
+    trackCharge[iTrack1] = digit[0]->TrackCharge(iTrack1);
+  }
+  nTracksTot = nTracks[0];
+  
+  // complement with contributing tracks on second cathod
+  for (Int_t iTrack2 = 0; iTrack2 < nTracks[1]; iTrack2++) {
+    Int_t trackId2 = digit[1]->Track(iTrack2);
+    // check if track exist
+    Bool_t trackExist = kFALSE;
+    for (Int_t iTrack1 = 0; iTrack1 < nTracks[0]; iTrack1++) {
+      if (trackId2 == trackId[iTrack1]) {
+       // complement existing track
+       trackCharge[iTrack1] += digit[1]->TrackCharge(iTrack2);
+       trackExist = kTRUE;
+       break;
+      }
+    }
+    // add the new track
+    if (!trackExist) {
+      trackId[nTracksTot] = trackId2;
+      trackCharge[nTracksTot] = digit[1]->TrackCharge(iTrack2);
+      nTracksTot++;
+    }
+  }
+  
+  // --- Find the most contributing track ---
+  Int_t mainTrackId = -1;
+  Float_t maxCharge = 0.;
+  for (Int_t iTrack = 0; iTrack < nTracksTot; iTrack++) {
+    if (trackCharge[iTrack] > maxCharge) {
+      mainTrackId = trackId[iTrack];
+      maxCharge = trackCharge[iTrack];
+    }
+  }
+  
+  delete[] trackId;
+  delete[] trackCharge;
+  
+  return mainTrackId;
+}
+
+//_____________________________________________________________________________
 void 
 AliMUONSimpleClusterServer::Print(Option_t*) const
 {
@@ -389,7 +464,7 @@ AliMUONSimpleClusterServer::Print(Option_t*) const
       {
         cout << Form("  -- Cathode %1d",cathode) << endl;
         
-        TClonesArray* padArray = PadArray(detElemId,cathode);
+        TObjArray* padArray = PadArray(detElemId,cathode);
         
         if (!padArray)
         {