]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MFT/AliMFTTrackerMU.cxx
Coordinates of SPD vertex
[u/mrichter/AliRoot.git] / MFT / AliMFTTrackerMU.cxx
index 75549ec6c6cb29d67a6af851a2b20ebedaf8e982..42338ca6c226269961bfd52fb1642dfa2b2cd788 100644 (file)
@@ -1,4 +1,4 @@
-/**************************************************************************
+/*************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 #include "AliMFTTrackerMU.h"
 #include "TMath.h"
 #include "AliRun.h"
+#include "AliRunLoader.h"
+#include "AliHeader.h"
+#include "AliGenEventHeader.h"
+#include "TArrayF.h"
 #include "AliMFT.h"
 #include "AliMUONTrackExtrap.h"
 #include "AliMUONTrack.h"
 #include "AliMUONESDInterface.h"
 #include "AliMuonForwardTrack.h"
 #include "AliMUONConstants.h"
+#include "TGrid.h"
 
 ClassImp(AliMFTTrackerMU)
 
@@ -61,9 +66,11 @@ AliMFTTrackerMU::AliMFTTrackerMU() :
   fMUONTrack(0),
   fCurrentTrack(0),
   fFinalBestCandidate(0),
-  fVertexErrorX(0.015),
-  fVertexErrorY(0.015),
-  fVertexErrorZ(0.010),
+  fXExtrapVertex(0),
+  fYExtrapVertex(0),
+  fZExtrapVertex(0),
+  fXExtrapVertexError(0),
+  fYExtrapVertexError(0),
   fBransonCorrection(kFALSE)
 {
 
@@ -126,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;
@@ -136,9 +147,9 @@ Int_t AliMFTTrackerMU::LoadClusters(TTree *cTree) {
 
 void AliMFTTrackerMU::UnloadClusters() {
   
-//   //--------------------------------------------------------------------
-//   // This function unloads MFT clusters
-//   //--------------------------------------------------------------------
+  //--------------------------------------------------------------------
+  // This function unloads MFT clusters
+  //--------------------------------------------------------------------
 
   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
     fMFTClusterArray[iPlane]      -> Clear("C");
@@ -168,15 +179,26 @@ Int_t AliMFTTrackerMU::Clusters2Tracks(AliESDEvent *event) {
 
   fESD = event;
 
-  // get the ITS primary vertex
+  fXExtrapVertex = 0;
+  fYExtrapVertex = 0;
+  fZExtrapVertex = 0;
+  fXExtrapVertexError = AliMFTConstants::fXVertexTolerance;
+  fYExtrapVertexError = AliMFTConstants::fYVertexTolerance;
 
-  Double_t vertex[3] = {0., 0., 0.};
+  // Imposing a fixed primary vertex, in order to keep memory of its position. Taking the primary vertex would imply the risk
+  // of loosing memory of its position when passing from ESD to AOD, due to possible refitting
   const AliESDVertex* esdVert = fESD->GetVertex(); 
   if (esdVert->GetNContributors() > 0 || !strcmp(esdVert->GetTitle(),"vertexer: smearMC")) {
-    esdVert->GetXYZ(vertex);
-    AliDebug(1,Form("found vertex (%e,%e,%e)",vertex[0],vertex[1],vertex[2]));
+    fXExtrapVertex = esdVert->GetX();
+    fYExtrapVertex = esdVert->GetY();
+    fZExtrapVertex = esdVert->GetZ();
+    fXExtrapVertexError = TMath::Max(AliMFTConstants::fXVertexTolerance, esdVert->GetXRes());
+    fYExtrapVertexError = TMath::Max(AliMFTConstants::fYVertexTolerance, esdVert->GetYRes());
+    AliInfo(Form("Found ESD vertex from %d contributors (%f +/- %f,  %f +/- %f,  %f)", 
+                esdVert->GetNContributors(),fXExtrapVertex,fXExtrapVertexError,fYExtrapVertex,fYExtrapVertexError,fZExtrapVertex));
   }
-  
+  else GetVertexFromMC();    
+
   //----------- Read ESD MUON tracks -------------------
 
   Int_t nTracksMUON = event->GetNumberOfMuonTracks();
@@ -200,7 +222,13 @@ Int_t AliMFTTrackerMU::Clusters2Tracks(AliESDEvent *event) {
     if (fMUONTrack) delete fMUONTrack;
     fMUONTrack = new AliMUONTrack();
 
-    AliMUONESDInterface::ESDToMUON(*esdTrack, *fMUONTrack, kTRUE);     // last arguments should be kTRUE or kFALSE??
+    AliMUONESDInterface::ESDToMUON(*esdTrack, *fMUONTrack, kFALSE);
+
+    if (!fMUONTrack->GetTrackParamAtCluster()->First()) {
+      AliInfo("Skipping track, no parameters available!!!");
+      iTrack++;
+      continue;
+    }
 
     // the track we are going to build, starting from fMUONTrack and adding the MFT clusters
     AliMuonForwardTrack *track = new ((*fCandidateTracks)[0]) AliMuonForwardTrack();
@@ -228,6 +256,8 @@ Int_t AliMFTTrackerMU::Clusters2Tracks(AliESDEvent *event) {
       
       Int_t nCandidates = fCandidateTracks->GetEntriesFast();
       for (Int_t iCandidate=0; iCandidate<nCandidates; iCandidate++) {
+
+       if (!(fCandidateTracks->UncheckedAt(iCandidate))) continue;
        fCurrentTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iCandidate);
 
        // if the old track is compatible with the new cluster, the track is updated and inserted as new track in the array 
@@ -270,6 +300,7 @@ Int_t AliMFTTrackerMU::Clusters2Tracks(AliESDEvent *event) {
 
     for (Int_t iFinalCandidate=0; iFinalCandidate<nFinalTracks; iFinalCandidate++) {
       
+      if (!(fCandidateTracks->UncheckedAt(iFinalCandidate))) continue;
       AliMuonForwardTrack *finalTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iFinalCandidate);
       Int_t nMFTClusters  = finalTrack->GetNMFTClusters();
 
@@ -298,26 +329,36 @@ Int_t AliMFTTrackerMU::Clusters2Tracks(AliESDEvent *event) {
 
       //----------------------- Save the information to the AliESDMuonGlobalTrack object
 
-      newTrack -> EvalKinem(vertex[2]);     // we evaluate the kinematics at the primary vertex
+      newTrack -> EvalKinem(AliMFTConstants::fZEvalKinem);
 
       AliDebug(2,"Creating a new Muon Global Track");
       AliESDMuonGlobalTrack *myESDTrack = event->NewMuonGlobalTrack();
       myESDTrack -> SetPxPyPz(newTrack->Px(), newTrack->Py(), newTrack->Pz());
       
-//       waiting for the commit of the new version of AliESDMuonGlobalTrack...
-
-//       myESDTrack -> SetChi2OverNdf(newTrack->GetChi2OverNdf());
-//       myESDTrack -> SetCharge(newTrack->GetCharge());
-//       myESDTrack -> SetMatchTrigger(newTrack->GetMatchTrigger());
-//       myESDTrack -> SetFirstTrackingPoint(newTrack->GetMFTCluster(0)->GetX(), newTrack->GetMFTCluster(0)->GetY(), newTrack->GetMFTCluster(0)->GetZ());
-//       myESDTrack -> SetXYAtVertex(newTrack->GetOffsetX(vertex[0], vertex[2]), newTrack->GetOffsetX(vertex[1], vertex[2]));
-//       myESDTrack -> SetRAtAbsorberEnd(newTrack->GetRAtAbsorberEnd());
-//       myESDTrack -> SetChi2MatchTrigger(esdTrack->GetChi2MatchTrigger());
-//       myESDTrack -> SetMuonClusterMap(esdTrack->GetMuonClusterMap());
-//       myESDTrack -> SetHitsPatternInTrigCh(esdTrack->GetHitsPatternInTrigCh());
-//       myESDTrack -> SetHitsPatternInTrigChTrk(esdTrack->GetHitsPatternInTrigChTrk());
-//       myESDTrack -> Connected(esdTrack->IsConnected());
-
+      myESDTrack -> SetLabel(newTrack->GetMCLabel());
+      myESDTrack -> SetChi2OverNdf(newTrack->GetChi2OverNdf());
+      myESDTrack -> SetCharge(newTrack->GetCharge());
+      myESDTrack -> SetMatchTrigger(newTrack->GetMatchTrigger());
+      myESDTrack -> SetNMFTClusters(newTrack->GetNMFTClusters());
+      myESDTrack -> SetNWrongMFTClustersMC(newTrack->GetNWrongClustersMC());
+      myESDTrack -> SetFirstTrackingPoint(newTrack->GetMFTCluster(0)->GetX(), newTrack->GetMFTCluster(0)->GetY(), newTrack->GetMFTCluster(0)->GetZ());
+      myESDTrack -> SetXYAtVertex(newTrack->GetOffsetX(0., AliMFTConstants::fZEvalKinem), newTrack->GetOffsetY(0., AliMFTConstants::fZEvalKinem));
+      myESDTrack -> SetRAtAbsorberEnd(newTrack->GetRAtAbsorberEnd());
+      myESDTrack -> SetCovariances(newTrack->GetTrackParamAtMFTCluster(0)->GetCovariances());
+      myESDTrack -> SetChi2MatchTrigger(esdTrack->GetChi2MatchTrigger());
+      myESDTrack -> SetMuonClusterMap(esdTrack->GetMuonClusterMap());
+      myESDTrack -> SetHitsPatternInTrigCh(esdTrack->GetHitsPatternInTrigCh());
+      myESDTrack -> SetHitsPatternInTrigChTrk(esdTrack->GetHitsPatternInTrigChTrk());
+      myESDTrack -> Connected(esdTrack->IsConnected());
+
+      ULong_t mftClusterPattern = 0;
+      for (Int_t iCluster=0; iCluster<newTrack->GetNMFTClusters(); iCluster++) {
+       AliMFTCluster *localCluster = newTrack->GetMFTCluster(iCluster);
+       mftClusterPattern |= 1 << localCluster->GetPlane();
+       mftClusterPattern |= IsCorrectMatch(localCluster, newTrack->GetMCLabel()) << (fNMaxPlanes + localCluster->GetPlane());
+      }
+      myESDTrack -> SetMFTClusterPattern(mftClusterPattern);
+      
       //---------------------------------------------------------------------------------
 
     }
@@ -379,19 +420,16 @@ Int_t AliMFTTrackerMU::FindClusterInPlane(Int_t planeId) {
     currentParamBack  = (*((AliMUONTrackParam*)(fMUONTrack->GetTrackParamAtCluster()->First())));
     currentParamForResearchFront = currentParamFront;
     currentParamForResearchBack  = currentParamBack;
-    Double_t xExtrap = gRandom->Gaus(0,fVertexErrorX);
-    Double_t yExtrap = gRandom->Gaus(0,fVertexErrorY);
-    Double_t zExtrap = gRandom->Gaus(0,fVertexErrorZ);
     if (fBransonCorrection) {
-      AliMUONTrackExtrap::ExtrapToVertex(&currentParamFront, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY); 
-      AliMUONTrackExtrap::ExtrapToVertex(&currentParamBack,  xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY); 
+      AliMUONTrackExtrap::ExtrapToVertex(&currentParamFront, fXExtrapVertex, fYExtrapVertex, fZExtrapVertex, fXExtrapVertexError, fYExtrapVertexError); 
+      AliMUONTrackExtrap::ExtrapToVertex(&currentParamBack,  fXExtrapVertex, fYExtrapVertex, fZExtrapVertex, fXExtrapVertexError, fYExtrapVertexError); 
     }
     else {
-      AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamFront, zExtrap);
-      AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamBack,  zExtrap);
+      AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamFront, fZExtrapVertex);
+      AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamBack,  fZExtrapVertex);
     }
-    AliMUONTrackExtrap::ExtrapToVertex(&currentParamForResearchFront, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY); 
-    AliMUONTrackExtrap::ExtrapToVertex(&currentParamForResearchBack,  xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY); 
+    AliMUONTrackExtrap::ExtrapToVertex(&currentParamForResearchFront, fXExtrapVertex, fYExtrapVertex, fZExtrapVertex, fXExtrapVertexError, fYExtrapVertexError); 
+    AliMUONTrackExtrap::ExtrapToVertex(&currentParamForResearchBack,  fXExtrapVertex, fYExtrapVertex, fZExtrapVertex, fXExtrapVertexError, fYExtrapVertexError); 
   }
   else {          // MFT planes others than the last one: mult. scattering correction because of the upstream MFT planes is performed
     currentParamFront = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
@@ -522,3 +560,150 @@ Double_t AliMFTTrackerMU::TryOneCluster(const AliMUONTrackParam &trackParam, Ali
 
 //=========================================================================================================================================
 
+Bool_t AliMFTTrackerMU::IsCorrectMatch(AliMFTCluster *cluster, Int_t labelMC) {
+
+  Bool_t result = kFALSE;
+
+  // check if the cluster belongs to the correct MC track
+
+  for (Int_t iTrack=0; iTrack<cluster->GetNMCTracks(); iTrack++) {
+    if (cluster->GetMCLabel(iTrack)==labelMC) {
+      result = kTRUE;
+      break;
+    }
+  }
+
+  return result;
+
+}
+
+//======================================================================================================================================
+
+void AliMFTTrackerMU::GetVertexFromMC() {
+
+  AliRunLoader *runLoader = AliRunLoader::Open("galice.root");
+  if (!runLoader) {
+    AliError("no run loader found in file galice.root");
+    return;
+  }
+
+  runLoader->CdGAFile();
+  runLoader->LoadgAlice();
+  runLoader->LoadHeader();
+  runLoader->GetEvent(gAlice->GetEvNumber());
+  
+  TArrayF vtx(3);
+  runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vtx);
+  AliInfo(Form("Primary vertex from MC found in (%f, %f, %f)\n",vtx[0], vtx[1], vtx[2]));
+
+  fXExtrapVertex = gRandom->Gaus(vtx[0], AliMFTConstants::fPrimaryVertexResX);
+  fYExtrapVertex = gRandom->Gaus(vtx[1], AliMFTConstants::fPrimaryVertexResY);
+  fZExtrapVertex = gRandom->Gaus(vtx[2], AliMFTConstants::fPrimaryVertexResZ);
+  fXExtrapVertexError = AliMFTConstants::fXVertexTolerance;
+  fYExtrapVertexError = AliMFTConstants::fYVertexTolerance;
+  AliInfo(Form("Set ESD vertex from MC (%f +/- %f,  %f +/- %f,  %f)", 
+              fXExtrapVertex,fXExtrapVertexError,fYExtrapVertex,fYExtrapVertexError,fZExtrapVertex));
+  
+}
+
+//======================================================================================================================================
+
+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];
+
+  }
+
+}
+
+//======================================================================================================================================
+