From 7bafb008a3b0b3ceb00933a81eb4f83348e80459 Mon Sep 17 00:00:00 2001
From: auras <antonio.uras@cern.ch>
Date: Thu, 24 Jul 2014 01:04:02 +0200
Subject: [PATCH] Modifications to MFT classes to manage embedding

---
 MFT/AliMFT.cxx              |  25 +++++++--
 MFT/AliMFT.h                |  18 ++++++-
 MFT/AliMFTAnalysisTools.cxx |  22 +++++---
 MFT/AliMFTConstants.h       |   4 +-
 MFT/AliMFTTrackerMU.cxx     | 104 ++++++++++++++++++++++++++++++++++++
 MFT/AliMFTTrackerMU.h       |   2 +
 6 files changed, 162 insertions(+), 13 deletions(-)

diff --git a/MFT/AliMFT.cxx b/MFT/AliMFT.cxx
index 0b1e35a6fd7..38e295abbc9 100644
--- a/MFT/AliMFT.cxx
+++ b/MFT/AliMFT.cxx
@@ -63,11 +63,18 @@ AliMFT::AliMFT():
   fChargeDispersion(25.e-4),
   fSingleStepForChargeDispersion(0),
   fNStepForChargeDispersion(4),
-  fDensitySupportOverSi(0.15)
+  fDensitySupportOverSi(0.15),
+  fFileNameForUnderyingEvent(0), 
+  fFileNameForPileUpEvents(0),
+  fNPileUpEvents(0), 
+  fUnderlyingEventID(-1)
 {
 
   // default constructor
 
+  for (Int_t iPileUp=0; iPileUp<AliMFTConstants::fNMaxPileUpEvents; iPileUp++) fPileUpEventsIDs[iPileUp] = -1;
+
+
 }
 
 //====================================================================================================================================================
@@ -86,9 +93,15 @@ AliMFT::AliMFT(const Char_t *name, const Char_t *title):
   fChargeDispersion(25.e-4),
   fSingleStepForChargeDispersion(0),
   fNStepForChargeDispersion(4),
-  fDensitySupportOverSi(0.15)
+  fDensitySupportOverSi(0.15),
+  fFileNameForUnderyingEvent(0), 
+  fFileNameForPileUpEvents(0),
+  fNPileUpEvents(0), 
+  fUnderlyingEventID(-1)
 {
 
+  for (Int_t iPileUp=0; iPileUp<AliMFTConstants::fNMaxPileUpEvents; iPileUp++) fPileUpEventsIDs[iPileUp] = -1;
+
   fNameGeomFile = "AliMFTGeometry.root";
 
   SetGeometry();
@@ -113,9 +126,15 @@ AliMFT::AliMFT(const Char_t *name, const Char_t *title, Char_t *nameGeomFile):
   fChargeDispersion(25.e-4),
   fSingleStepForChargeDispersion(0),
   fNStepForChargeDispersion(4),
-  fDensitySupportOverSi(0.15)
+  fDensitySupportOverSi(0.15),
+  fFileNameForUnderyingEvent(0), 
+  fFileNameForPileUpEvents(0),
+  fNPileUpEvents(0), 
+  fUnderlyingEventID(-1)
 {
 
+  for (Int_t iPileUp=0; iPileUp<AliMFTConstants::fNMaxPileUpEvents; iPileUp++) fPileUpEventsIDs[iPileUp] = -1;
+
   fNameGeomFile = nameGeomFile;
 
   SetGeometry();
diff --git a/MFT/AliMFT.h b/MFT/AliMFT.h
index fb3e50f7c75..ae612135b0b 100644
--- a/MFT/AliMFT.h
+++ b/MFT/AliMFT.h
@@ -102,6 +102,19 @@ public:
 
   void SetDensitySupportOverSi(Double_t density) { if (density>1e-6) fDensitySupportOverSi=density; else fDensitySupportOverSi=1e-6; }
 
+  //--------- for underlying and pile-up events --------------------
+
+  void SetFileNameForUnderlyingEvent(TString fileName) { fFileNameForUnderyingEvent += fileName; }
+  void SetFileNameForPileUpEvents(TString fileName)    { fFileNameForPileUpEvents   += fileName; }
+
+  void SetUnderlyingEventID(Short_t eventID) { fUnderlyingEventID = eventID; }
+  void SetPileUpEventID(Short_t i, Short_t eventID) { if (i>=0 && i<AliMFTConstants::fNMaxPileUpEvents) fPileUpEventsIDs[i] = eventID; }
+
+  const Char_t* GetFileNameForUnderlyingEvent() { return fFileNameForUnderyingEvent; }
+  const Char_t* GetFileNameForPileUpEvents()    { return fFileNameForPileUpEvents; }
+  Short_t GetUnderlyingEventID() { return fUnderlyingEventID; }
+  Short_t GetPileUpEventID(Short_t i) { if (i>=0 && i<AliMFTConstants::fNMaxPileUpEvents) return fPileUpEventsIDs[i]; else return -1; }
+  
 protected:
 
   static const Int_t fNMaxPlanes = AliMFTConstants::fNMaxPlanes;        // max number of MFT planes
@@ -125,13 +138,16 @@ protected:
   Int_t fNStepForChargeDispersion;
 
   Double_t fDensitySupportOverSi;
+
+  TString fFileNameForUnderyingEvent, fFileNameForPileUpEvents;
+  Short_t fNPileUpEvents, fUnderlyingEventID, fPileUpEventsIDs[AliMFTConstants::fNMaxPileUpEvents];
   
 private:
 
   AliMFT (const AliMFT& mft);             // dummy copy constructor
   AliMFT &operator=(const AliMFT& mft);   // dummy assignment operator
   
-  ClassDef(AliMFT,1) 
+  ClassDef(AliMFT,2) 
     
 };
 
diff --git a/MFT/AliMFTAnalysisTools.cxx b/MFT/AliMFTAnalysisTools.cxx
index abce117eba0..34fef4b7edd 100644
--- a/MFT/AliMFTAnalysisTools.cxx
+++ b/MFT/AliMFTAnalysisTools.cxx
@@ -31,6 +31,7 @@
 #include "TMath.h"
 #include "AliLog.h"
 #include "TObjArray.h"
+#include "TDecompLU.h"
 
 #include "AliMFTAnalysisTools.h"
 
@@ -152,16 +153,21 @@ Double_t AliMFTAnalysisTools::GetAODMuonWeightedOffset(AliAODTrack *muon, Double
   covCoordinates(1,0) = cov(2,0);
   covCoordinates(1,1) = cov(2,2);
   
-  TMatrixD covCoordinatesInverse = covCoordinates.Invert();
+  if (TDecompLU::InvertLU(covCoordinates,covCoordinates.GetTol(),0)) {
 
-  Double_t dX = xy[0] - xv;
-  Double_t dY = xy[1] - yv;
-  
-  Double_t weightedOffset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) + 
-					     dY*dY*covCoordinatesInverse(1,1) + 
-					     2.*dX*dY*covCoordinatesInverse(0,1)));
+    TMatrixD covCoordinatesInverse = covCoordinates;
+    Double_t dX = xy[0] - xv;
+    Double_t dY = xy[1] - yv;
+    
+    Double_t weightedOffset = TMath::Sqrt(0.5*(dX*dX*covCoordinatesInverse(0,0) +
+					       dY*dY*covCoordinatesInverse(1,1) +
+					       2.*dX*dY*covCoordinatesInverse(0,1)));
+    
+    return weightedOffset;
 
-  return weightedOffset;
+  }
+  
+  return -999;
 
 }
 
diff --git a/MFT/AliMFTConstants.h b/MFT/AliMFTConstants.h
index e96553d51f0..1cbb3b4cffe 100644
--- a/MFT/AliMFTConstants.h
+++ b/MFT/AliMFTConstants.h
@@ -58,13 +58,15 @@ public:
   static const Double_t fMisalignmentMagnitude;   // Expected misalignment magnitude (for MC, waiting for OCDB)
 
   static const Int_t fNMaxMuonsForPCA = 10;
+  static const Int_t fNMaxPileUpEvents = 5;         // Max events to consider for pile-up studies
+  static const Int_t fLabelOffsetMC = 10000000;     // Offset to be added to MC labels of tracks from merged underlying and pile-up events
 
 protected:
 
   AliMFTConstants() : TObject() {}
   virtual ~AliMFTConstants(){}
 
-  ClassDef(AliMFTConstants, 3)    // MFT global constants 
+  ClassDef(AliMFTConstants, 4)    // MFT global constants 
 
 };
 	
diff --git a/MFT/AliMFTTrackerMU.cxx b/MFT/AliMFTTrackerMU.cxx
index 422bb9c1a43..0204c6eab8c 100644
--- a/MFT/AliMFTTrackerMU.cxx
+++ b/MFT/AliMFTTrackerMU.cxx
@@ -43,6 +43,7 @@
 #include "AliMUONESDInterface.h"
 #include "AliMuonForwardTrack.h"
 #include "AliMUONConstants.h"
+#include "TGrid.h"
 
 ClassImp(AliMFTTrackerMU)
 
@@ -132,6 +133,10 @@ Int_t AliMFTTrackerMU::LoadClusters(TTree *cTree) {
   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
     AliInfo(Form("plane %02d: nClusters = %d\n", iPlane, fMFTClusterArray[iPlane]->GetEntries()));
   }
+
+  AddClustersFromUnderlyingEvent();
+  AddClustersFromPileUpEvents();
+
   SeparateFrontBackClusters();
 
   return 0;
@@ -599,3 +604,102 @@ void AliMFTTrackerMU::GetVertexFromMC() {
 
 //======================================================================================================================================
 
+void AliMFTTrackerMU::AddClustersFromUnderlyingEvent() {
+
+  AliInfo("Adding clusters from underlying event");
+
+  if (!fMFT) return;
+
+  TGrid::Connect("alien://");
+
+  TFile* fileWithClustersToAdd = TFile::Open(fMFT->GetFileNameForUnderlyingEvent());
+  if (!fileWithClustersToAdd) return;
+  if (!(fileWithClustersToAdd->IsOpen())) return;
+  if (!(fileWithClustersToAdd->cd(Form("Event%d",fMFT->GetUnderlyingEventID())))) return;
+
+  TClonesArray *recPointsPerPlaneToAdd[AliMFTConstants::fNMaxPlanes] = {0};
+  TTree *treeIn = 0;
+
+  for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) recPointsPerPlaneToAdd[iPlane] = new TClonesArray("AliMFTCluster");
+
+  treeIn = (TTree*) gDirectory->Get("TreeR");
+
+  for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
+    if (!(treeIn->GetBranch(Form("Plane_%02d",iPlane)))) {
+      for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) delete recPointsPerPlaneToAdd[jPlane];
+      return;
+    }
+    else treeIn->SetBranchAddress(Form("Plane_%02d",iPlane), &(recPointsPerPlaneToAdd[iPlane]));
+  }
+
+  treeIn -> GetEntry(0);
+
+  for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
+    printf("plane %d -> before = %d ",iPlane,fMFTClusterArray[iPlane]->GetEntries());
+    Int_t nClusters = recPointsPerPlaneToAdd[iPlane]->GetEntries();
+    for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
+      AliMFTCluster *newCluster = (AliMFTCluster*) recPointsPerPlaneToAdd[iPlane]->At(iCluster);
+      for (Int_t iTrack=0; iTrack<newCluster->GetNMCTracks(); iTrack++) newCluster->SetMCLabel(iTrack, newCluster->GetMCLabel(iTrack)+AliMFTConstants::fLabelOffsetMC);
+      new ((*fMFTClusterArray[iPlane])[fMFTClusterArray[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
+    }
+    printf("after = %d\n",fMFTClusterArray[iPlane]->GetEntries());
+  }
+
+  for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) delete recPointsPerPlaneToAdd[jPlane];
+
+}
+
+//======================================================================================================================================
+
+void AliMFTTrackerMU::AddClustersFromPileUpEvents() {
+
+  AliInfo("Adding clusters from pile-up event(s)");
+
+  if (!fMFT) return;
+
+  TGrid::Connect("alien://");
+
+  TFile* fileWithClustersToAdd = TFile::Open(fMFT->GetFileNameForPileUpEvents());
+  if (!fileWithClustersToAdd) return;
+  if (!(fileWithClustersToAdd->IsOpen())) return;
+
+  TClonesArray *recPointsPerPlaneToAdd[AliMFTConstants::fNMaxPlanes] = {0};
+  TTree *treeIn = 0;
+
+  for (Int_t iPileUp=0; iPileUp<AliMFTConstants::fNMaxPileUpEvents; iPileUp++) {
+    
+    if (!(fileWithClustersToAdd->cd(Form("Event%d",fMFT->GetPileUpEventID(iPileUp))))) continue;
+    
+    for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) recPointsPerPlaneToAdd[iPlane] = new TClonesArray("AliMFTCluster");
+    
+    treeIn = (TTree*) gDirectory->Get("TreeR");
+    
+    for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
+      if (!(treeIn->GetBranch(Form("Plane_%02d",iPlane)))) {
+	for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) delete recPointsPerPlaneToAdd[jPlane];
+	return;
+      }
+      else treeIn->SetBranchAddress(Form("Plane_%02d",iPlane), &(recPointsPerPlaneToAdd[iPlane]));
+    }
+
+    treeIn -> GetEntry(0);
+    
+    for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
+      printf("plane %d -> before = %d ",iPlane,fMFTClusterArray[iPlane]->GetEntries());
+      Int_t nClusters = recPointsPerPlaneToAdd[iPlane]->GetEntries();
+      for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
+	AliMFTCluster *newCluster = (AliMFTCluster*) recPointsPerPlaneToAdd[iPlane]->At(iCluster);
+	for (Int_t iTrack=0; iTrack<newCluster->GetNMCTracks(); iTrack++) newCluster->SetMCLabel(iTrack, newCluster->GetMCLabel(iTrack)+AliMFTConstants::fLabelOffsetMC);
+	new ((*fMFTClusterArray[iPlane])[fMFTClusterArray[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
+      }
+      printf("after = %d\n",fMFTClusterArray[iPlane]->GetEntries());
+    }
+
+    for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) delete recPointsPerPlaneToAdd[jPlane];
+
+  }
+
+}
+
+//======================================================================================================================================
+
diff --git a/MFT/AliMFTTrackerMU.h b/MFT/AliMFTTrackerMU.h
index f4c0ddc2ebd..310852a7ad6 100644
--- a/MFT/AliMFTTrackerMU.h
+++ b/MFT/AliMFTTrackerMU.h
@@ -64,6 +64,8 @@ public:
   /// Dummy implementation
   virtual AliCluster *GetCluster(Int_t /*index*/) const {return 0;}
 
+  void AddClustersFromUnderlyingEvent();
+  void AddClustersFromPileUpEvents();
 
 protected:
 
-- 
2.43.0