fix finding of pad neighbours; remove methods to write them in OCDB
[u/mrichter/AliRoot.git] / MUON / AliMUONSimpleClusterServer.cxx
index c0e34df..b6bfe99 100644 (file)
 #include "AliMpDEIterator.h"
 #include "AliMpDEManager.h"
 #include "AliMpExMap.h"
+#include "AliMpExMapIterator.h"
 #include "AliMpPad.h"
 #include "AliMpSegmentation.h"
 #include "AliMpVSegmentation.h"
 #include <Riostream.h>
-#include <TClonesArray.h>
+#include <TObjArray.h>
 #include <TString.h>
 #include <float.h>
 
@@ -48,6 +49,8 @@
 /// 
 /// \author Laurent Aphecetche, Subatech
 
+using std::endl;
+using std::cout;
 /// \cond CLASSIMP  
 ClassImp(AliMUONSimpleClusterServer)
 /// \endcond
@@ -57,17 +60,18 @@ 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,
                                                        const AliMUONGeometryTransformer& transformer)
-: AliMUONVClusterServer(), 
+: AliMUONVClusterServer(),
+  fDigitStore(0x0),
   fClusterFinder(clusterFinder),
   fkTransformer(transformer),
   fPads(),
@@ -80,6 +84,8 @@ AliMUONSimpleClusterServer::AliMUONSimpleClusterServer(AliMUONVClusterFinder* cl
     fPads[0] = new AliMpExMap;
     fPads[1] = new AliMpExMap;
     
+    fPadsIterator[0] = fPads[0]->CreateIterator();
+    fPadsIterator[1] = fPads[1]->CreateIterator();
 }
 
 //_____________________________________________________________________________
@@ -89,6 +95,8 @@ AliMUONSimpleClusterServer::~AliMUONSimpleClusterServer()
   delete fClusterFinder;
   delete fPads[0];
   delete fPads[1];
+  delete fPadsIterator[0];
+  delete fPadsIterator[1];
   delete fBypass;
 }
 
@@ -105,7 +113,7 @@ 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));
+  AliCodeTimerAuto(Form("Chamber %d",chamberId),0);
   
   if ( fTriggerTrackStore && chamberId >= 6 ) 
   {
@@ -131,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 ) || 
@@ -150,12 +158,16 @@ 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->SetChargeHints(recoParam->LowestPadCharge(),
+                                       recoParam->LowestClusterCharge());
+        
+        if ( fClusterFinder->NeedSegmentation() )
+        {
           fClusterFinder->Prepare(detElemId,pads,deArea,seg);
         }
         else
@@ -203,10 +215,16 @@ AliMUONSimpleClusterServer::Clusterize(Int_t chamberId,
           rawCluster->SetXYZ(xg, yg, zg);
           rawCluster->SetErrXY(recoParam->GetDefaultNonBendingReso(chamberId),recoParam->GetDefaultBendingReso(chamberId));
           
-          AliDebug(1,Form("Adding RawCluster detElemId %4d mult %2d charge %e (xl,yl,zl)=(%e,%e,%e) (xg,yg,zg)=(%e,%e,%e)",
+          // 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) label %d",
                           detElemId,rawCluster->GetNDigits(),rawCluster->GetCharge(),
                           cluster->Position().X(),cluster->Position().Y(),0.0,
-                          xg,yg,zg));
+                          xg,yg,zg,rawCluster->GetMCLabel()));
         }
       }
     }
@@ -214,7 +232,7 @@ AliMUONSimpleClusterServer::Clusterize(Int_t chamberId,
   }
   
   AliDebug(1,Form("chamberId = %2d NofClusters after = %d",chamberId,fNCluster));
-
+  
   return nofAddedClusters;
 }
 
@@ -228,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);
   
   fkTransformer.Global2Local(detElemId,
-                             globalArea.Position().X(),globalArea.Position().Y(),zg,
+                             globalArea.GetPositionX(),globalArea.GetPositionY(),zg,
                              xl,yl,zl);
   
-  localArea = AliMpArea(TVector2(xl,yl), globalArea.Dimensions());
+  localArea = AliMpArea(xl,yl, globalArea.GetDimensionX(), globalArea.GetDimensionY());
 }
 
 //_____________________________________________________________________________
@@ -278,12 +301,12 @@ 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));
 }
 
 //_____________________________________________________________________________
@@ -299,43 +322,124 @@ AliMUONSimpleClusterServer::UseTriggerTrackStore(AliMUONVTriggerTrackStore* trac
 
 //_____________________________________________________________________________
 void 
-AliMUONSimpleClusterServer::UseDigits(TIter& next)
+AliMUONSimpleClusterServer::UseDigits(TIter& next, AliMUONVDigitStore* digitStore)
 {
   /// Convert digitStore into two arrays of AliMUONPads
-
-  fPads[0]->Clear();
-  fPads[1]->Clear();
   
+  fDigitStore = digitStore;
+  
+  // 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
 {
@@ -360,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)
         {