Services added for the MFT (R. Tieulent). Class for the Muon Global Tracking added.
authorAntonio Uras <antonio.uras@cern.ch>
Tue, 28 Jan 2014 17:29:38 +0000 (18:29 +0100)
committerAntonio Uras <antonio.uras@cern.ch>
Tue, 28 Jan 2014 17:29:38 +0000 (18:29 +0100)
MFT/AliMFT.cxx
MFT/AliMFT.h
MFT/AliMFTReconstructor.cxx
MFT/AliMFTReconstructor.h
MFT/AliMFTTrackerMU.cxx [new file with mode: 0644]
MFT/AliMFTTrackerMU.h [new file with mode: 0644]
MFT/AliMuonForwardTrackFinder.cxx
MFT/AliMuonForwardTrackFinder.h
MFT/CMakelibMFTrec.pkg
MFT/MFTrecLinkDef.h

index f20a504..4e46ca4 100644 (file)
@@ -142,9 +142,23 @@ void AliMFT::CreateMaterials() {
 
   AliInfo("Start MFT materials");
 
-  // data from PDG booklet 2002     density [gr/cm^3] rad len [cm] abs len [cm]    
-  Float_t   aAir[4]={12,14,16,36} ,   zAir[4]={6,7,8,18} ,   wAir[4]={0.000124,0.755267,0.231781,0.012827} , dAir=0.00120479; Int_t nAir=4;   // Air mixture
-  Float_t   aSi = 28.085 ,            zSi   = 14 ,           dSi   =  2.329    ,   radSi   =  21.82/dSi , absSi   = 108.4/dSi  ;              // Silicon
+  //---------------------- Materials and Mixtures ----------------------------------------------------------------------------------------------------
+
+  // data from PDG booklet 2002            density [gr/cm^3]    rad len [cm]             abs len [cm]    
+                                                                                        
+  Float_t   aSi = 28.085 ,  zSi   = 14. ,  dSi   = 2.329 ,      radSi   =  21.82/dSi ,   absSi   = 108.4/dSi  ;         // Silicon
+  Float_t   aCarb = 12.01 , zCarb =  6. ,  dCarb = 2.265 ,      radCarb =  18.8 ,        absCarb = 49.9  ;              // Carbon
+  Float_t   aAlu = 26.98 ,  zAlu  = 13. ,  dAlu  = 2.70  ,      radAlu  =  8.897 ,       absAlu  = 39.70  ;             // Aluminum
+
+  const Int_t nAir   = 4;
+  const Int_t nWater = 2;
+  const Int_t nSiO2  = 2;
+
+  Float_t   aAir[nAir]     = {12,14,16,40} ,      zAir[nAir]     = {6,7,8,18} ,   wAir[nAir]     = {0.000124,0.755267,0.231781,0.012827} , dAir=0.00120479; // Air mixture
+  Float_t   aWater[nWater] = {1.00794,15.9994} ,  zWater[nWater] = {1,8} ,        wWater[nWater] = {0.111894,0.888106} ,                   dWater=1.;       // Water mixture
+  Float_t   aSiO2[nSiO2]   = {15.9994,28.0855} ,  zSiO2[nSiO2]   = {8.,14.} ,     wSiO2[nSiO2]   = {0.532565,0.467435} ,                   dSiO2=2.20;      // SiO2 mixture
+
+  //---------------------------------------------------------------------------------------------------------------------------------------------------
 
   Int_t   matId  = 0;                        // tmp material id number
   Int_t   unsens = 0, sens=1;                // sensitive or unsensitive medium
@@ -177,7 +191,19 @@ void AliMFT::CreateMaterials() {
 
   AliMaterial(++matId, "Support", aSi,   zSi,    dSi*fDensitySupportOverSi, radSi/fDensitySupportOverSi, absSi/fDensitySupportOverSi);  
   AliMedium(kSupport,  "Support", matId, unsens, isxfld,  sxmgmx,    tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
-    
+  
+  AliMaterial(++matId, "Carbon",   aCarb,   zCarb,    dCarb,    radCarb,  absCarb );
+  AliMedium(kCarbon,  "Carbon", matId, unsens, isxfld,  sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+  
+  AliMaterial(++matId, "Alu",   aAlu,   zAlu,    dAlu,    radAlu,  absAlu );
+  AliMedium(kAlu,  "Alu", matId, unsens, isxfld,  sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+  
+  AliMixture(++matId,"Water", aWater,  zWater,   dWater,   nWater,   wWater);
+  AliMedium(kWater,    "Water", matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
+  
+  AliMixture(++matId,"SiO2", aSiO2,  zSiO2,   dSiO2,   nSiO2,   wSiO2);
+  AliMedium(kSiO2,    "SiO2", matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
+
   AliInfo("End MFT materials");
           
 }
@@ -274,22 +300,10 @@ void AliMFT::StepManager() {
   hit.SetMomentum(momentum);
   hit.SetStatus(status);
   hit.SetEloss(gMC->Edep());
-  //  hit.SetShunt(GetIshunt());
-//   if (gMC->IsTrackEntering()) {
-//     hit.SetStartPosition(position);
-//     hit.SetStartTime(gMC->TrackTime());
-//     hit.SetStartStatus(status);
-//     return; // don't save entering hit.
-//   } 
 
   // Fill hit structure with this new hit.
   new ((*fHits)[fNhits++]) AliMFTHit(hit);
 
-  // Save old position... for next hit.
-//   hit.SetStartPosition(position);
-//   hit.SetStartTime(gMC->TrackTime());
-//   hit.SetStartStatus(status);
-
   return;
 
 }
@@ -306,9 +320,91 @@ TGeoVolumeAssembly* AliMFT::CreateVol() {
   TGeoMedium *silicon = gGeoManager->GetMedium("MFT_Si");
   TGeoMedium *readout = gGeoManager->GetMedium("MFT_Readout");
   TGeoMedium *support = gGeoManager->GetMedium("MFT_Support");
-  for (Int_t iPar=0; iPar<8; iPar++) AliInfo(Form("silicon->GetParam(%d) = %f", iPar, silicon->GetParam(iPar)));
-  for (Int_t iPar=0; iPar<8; iPar++) AliInfo(Form("readout->GetParam(%d) = %f", iPar, readout->GetParam(iPar)));
-  for (Int_t iPar=0; iPar<8; iPar++) AliInfo(Form("support->GetParam(%d) = %f", iPar, support->GetParam(iPar)));
+  TGeoMedium *carbon  = gGeoManager->GetMedium("MFT_Carbon");
+  TGeoMedium *alu     = gGeoManager->GetMedium("MFT_Alu");
+  TGeoMedium *water   = gGeoManager->GetMedium("MFT_Water");
+  TGeoMedium *si02    = gGeoManager->GetMedium("MFT_SiO2");
+
+  // ---- Cage & Services Description --------------------------------------------------------------------
+
+  // R. Tieulent - 17/01/2014 - Basic description for ITS/TPC matching studies
+
+  TGeoVolumeAssembly *cageNservices = new TGeoVolumeAssembly("MFT_cageNservices");
+  
+  // cage definition
+  Float_t cage_dz   = 150./2.;
+  Float_t cage_rMin = 50.;
+  Float_t cage_rMax = cage_rMin + 0.188; // 1% of X0 for Carbon
+  
+  TGeoVolume *cage = gGeoManager->MakeTube("MFT_cage", carbon, cage_rMin, cage_rMax, cage_dz);
+  cageNservices->AddNode(cage,1,new TGeoTranslation(0., 0., 0. ));
+
+  // Services definition
+  TGeoVolumeAssembly *services = new TGeoVolumeAssembly("MFT_services");
+
+  // Aluminum bus-Bar
+  Float_t busBar_dz = 150.;
+  Float_t busBar_thick = 0.1 ;
+  Float_t busBar_large = 1.;
+  
+  TGeoVolume *aluBusBar = gGeoManager->MakeBox("MFT_busBar", alu, busBar_large/2., busBar_thick/2., busBar_dz/2.);
+
+  Int_t nBusBar = 30;
+  Float_t dPhi_busBar = 2.*TMath::Pi() / nBusBar;
+  Float_t dShift = cage_rMin - busBar_thick/2.;
+  
+  TGeoRotation *rot = 0;
+
+  for (Int_t iBusBar=0; iBusBar<nBusBar; iBusBar++) {
+    Float_t phi = dPhi_busBar*iBusBar;
+    Float_t xp  = dShift*TMath::Cos(phi);
+    Float_t yp  = dShift*TMath::Sin(phi);
+    rot = new TGeoRotation();
+    rot -> RotateZ(phi*TMath::RadToDeg()+90.);
+    services -> AddNode(aluBusBar,iBusBar+1,new TGeoCombiTrans(xp,yp,0,rot));
+  }
+  
+  // Cooling Water Pipes
+  Float_t cooling_dz = 150.;
+  Float_t cooling_r = 0.3/2. ; // 3mm in diameter
+  
+  TGeoVolume *cooling = gGeoManager->MakeTube("MFT_cooling", water, 0., cooling_r, cooling_dz/2.);
+  
+  Int_t nCooling = 18;
+  dShift = cage_rMin - cooling_r ;
+  Float_t phi0 = 0.02;
+  
+  for (Int_t iCooling=0; iCooling<nCooling; iCooling++) {
+    Float_t phi;
+    if (iCooling < nCooling/2) phi = dPhi_busBar*(iCooling+3) + phi0;
+    else                       phi = dPhi_busBar*(iCooling+9) + phi0;
+    Float_t xp = dShift*TMath::Cos(phi);
+    Float_t yp = dShift*TMath::Sin(phi);
+    services -> AddNode(cooling,iCooling+1,new TGeoTranslation(xp,yp,0 ));
+  }
+  
+  // Optical Fibers
+  Float_t fiber_dz = 150.;
+  Float_t fiber_r = 0.0125/2. ; // 0.125mm in diameter
+  
+  TGeoVolume *fiber = gGeoManager->MakeTube("MFT_fiber", si02, 0., fiber_r, fiber_dz/2.);
+  
+  Int_t nFiber = 340;
+  dShift = cage_rMin - fiber_r - cooling_r;
+  phi0 = 0.03;
+  
+  for (Int_t iFiber=0; iFiber<nFiber; iFiber++) {
+    Float_t phi = dPhi_busBar*(Int_t)(iFiber/11+1) - phi0-(iFiber%11)*2.*TMath::ATan(fiber_r/dShift);
+    Float_t xp  = dShift*TMath::Cos(phi);
+    Float_t yp  = dShift*TMath::Sin(phi);
+    services->AddNode(fiber,iFiber+1,new TGeoTranslation(xp,yp,0 ));
+  }
+
+  cageNservices->AddNode(services,1,new TGeoTranslation(0., 0., 0. ));
+
+  vol->AddNode(cageNservices,1,new TGeoTranslation(0., 0., 0. ));
+  
+  // ---------------- Planes description ------------------------------------------------------------------
 
   Double_t origin[3] = {0};
 
index af6a618..1d97796 100644 (file)
@@ -82,7 +82,7 @@ public:
   
   AliMFTSegmentation* GetSegmentation() const { return fSegmentation; }
 
-  enum EMedia{kAir, kSi, kReadout, kSupport};  // media IDs used in CreateMaterials  
+  enum EMedia{kAir, kSi, kReadout, kSupport, kCarbon, kAlu, kWater, kSiO2};  // media IDs used in CreateMaterials
     
   // Geometry/segmentation creation part
   TGeoVolumeAssembly* CreateVol();
index 5d0716a..b07dfa0 100644 (file)
@@ -91,7 +91,7 @@ void AliMFTReconstructor::Init() {
     ((TClonesArray*)fDigits->At(iPlane))->SetOwner(kTRUE);
   }
   
-  AliInfo("    ************* Using the MFT reconstructor! ****** ");
+  AliInfo(" ************* Using the MFT reconstructor! ************ ");
   
   return;
 
@@ -159,3 +159,30 @@ void AliMFTReconstructor::Reconstruct(TTree *digitsTree, TTree *clustersTree) co
 }
 
 //====================================================================================================================================================
+
+AliTracker* AliMFTReconstructor::CreateTracker() const {
+
+  // Create a MFT tracker: global MUON+MFT tracks
+
+  AliMFTTrackerMU *tracker = new AliMFTTrackerMU();   // tracker for muon tracks (MFT + MUON) 
+
+  return tracker;
+  
+}
+
+//====================================================================================================================================================
+
+AliTracker* AliMFTReconstructor::CreateTrackleter() const {
+
+  AliInfo("Not implemented");
+
+  // Create a MFT tracker: standalone MFT tracks
+
+  //  AliMFTTrackerSA *tracker = new AliMFTTrackerSA();   // tracker for MFT tracks
+
+  return NULL;
+  
+}
+
+//====================================================================================================================================================
+
index 268ed3c..e12b3a0 100644 (file)
@@ -16,6 +16,9 @@
 #include "TTree.h"
 #include "AliMFTSegmentation.h"
 #include "AliReconstructor.h"
+#include "AliTracker.h"
+#include "AliVertexer.h"
+#include "AliMFTTrackerMU.h"
 #include "AliMFTClusterFinder.h"
 
 //====================================================================================================================================================
@@ -37,6 +40,9 @@ public:
   virtual void  Reconstruct(TTree *digitsTree, TTree *clustersTree) const; 
   virtual void  Reconstruct(AliRawReader* /*rawdata*/, TTree* /*clustersTree*/) const { AliInfo("Not implemented"); } 
 
+  virtual AliTracker* CreateTracker()  const;
+  virtual AliTracker* CreateTrackleter()  const;
+
   //  static const AliMFTRecoParam* GetRecoParam() { return dynamic_cast<const AliMFTRecoParam*>(AliReconstructor::GetRecoParam(0)); }
 
 private:
diff --git a/MFT/AliMFTTrackerMU.cxx b/MFT/AliMFTTrackerMU.cxx
new file mode 100644 (file)
index 0000000..ee85e07
--- /dev/null
@@ -0,0 +1,471 @@
+/**************************************************************************
+* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+*                                                                        *
+* Author: The ALICE Off-line Project.                                    *
+* Contributors are mentioned in the code where appropriate.              *
+*                                                                        *
+* Permission to use, copy, modify and distribute this software and its   *
+* documentation strictly for non-commercial purposes is hereby granted   *
+* without fee, provided that the above copyright notice appears in all   *
+* copies and that both the copyright notice and this permission notice   *
+* appear in the supporting documentation. The authors make no claims     *
+* about the suitability of this software for any purpose. It is          *
+* provided "as is" without express or implied warranty.                  *
+**************************************************************************/
+
+//====================================================================================================================================================
+//
+//                          MFT tracker
+//
+// Class for the creation of the "global muon tracks" built from the tracks reconstructed in the 
+// muon spectrometer and the clusters of the Muon Forward Tracker
+//
+// Contact author: antonio.uras@cern.ch
+//
+//====================================================================================================================================================
+
+#include "TTree.h"
+#include "AliLog.h"
+#include "AliGeomManager.h"
+#include "AliESDEvent.h"
+#include "AliESDMuonTrack.h"
+#include "AliESDMuonGlobalTrack.h"
+#include "AliMFTTrackerMU.h"
+#include "TMath.h"
+#include "AliRun.h"
+#include "AliMFT.h"
+#include "AliMUONTrackExtrap.h"
+#include "AliMUONTrack.h"
+#include "AliMUONESDInterface.h"
+#include "AliMuonForwardTrack.h"
+
+ClassImp(AliMFTTrackerMU)
+
+const Double_t AliMFTTrackerMU::fRadLengthSi = AliMFTConstants::fRadLengthSi;
+
+//====================================================================================================================================================
+
+AliMFTTrackerMU::AliMFTTrackerMU() : 
+  AliTracker(),
+  fESD(0),
+  fMFT(0),
+  fSegmentation(0),
+  fNPlanesMFT(0),
+  fNPlanesMFTAnalyzed(0),
+  fSigmaClusterCut(0),
+  fScaleSigmaClusterCut(1.),
+  fNMaxMissingMFTClusters(0),
+  fGlobalTrackingDiverged(kFALSE),
+  fCandidateTracks(0),
+  fMUONTrack(0),
+  fCurrentTrack(0),
+  fFinalBestCandidate(0),
+  fVertexErrorX(0.015),
+  fVertexErrorY(0.015),
+  fVertexErrorZ(0.010),
+  fBransonCorrection(kFALSE)
+{
+
+  //--------------------------------------------------------------------
+  // This is the AliMFTTrackerMU constructor
+  //--------------------------------------------------------------------
+
+  fMFT = (AliMFT*) gAlice->GetDetector("MFT");
+  fSegmentation = fMFT->GetSegmentation();
+  SetNPlanesMFT(fSegmentation->GetNPlanes());
+  AliMUONTrackExtrap::SetField();                 // set the magnetic field for track extrapolations
+
+  for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
+    fMFTClusterArray[iPlane]      = 0;
+    fMFTClusterArrayFront[iPlane] = new TClonesArray("AliMFTCluster");
+    fMFTClusterArrayBack[iPlane]  = new TClonesArray("AliMFTCluster");
+    fIsPlaneMandatory[iPlane]     = kFALSE;
+  }
+
+}
+
+//====================================================================================================================================================
+
+AliMFTTrackerMU::~AliMFTTrackerMU() {
+
+  // destructor
+
+  for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
+    delete fMFTClusterArray[iPlane];
+    delete fMFTClusterArrayFront[iPlane];
+    delete fMFTClusterArrayBack[iPlane];
+  }
+
+}
+
+//====================================================================================================================================================
+
+Int_t AliMFTTrackerMU::LoadClusters(TTree *cTree) {
+
+  //--------------------------------------------------------------------
+  // This function loads the MFT clusters
+  //--------------------------------------------------------------------
+  if (!cTree->GetEvent()) return kFALSE;
+  for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
+    AliDebug(1, Form("plane %02d: nClusters = %d\n", iPlane, (fMFT->GetRecPointsList(iPlane))->GetEntries()));
+    fMFTClusterArray[iPlane] = fMFT->GetRecPointsList(iPlane);
+  }
+  SeparateFrontBackClusters();
+
+  return 0;
+
+}
+
+//====================================================================================================================================================
+
+void AliMFTTrackerMU::UnloadClusters() {
+  
+  //--------------------------------------------------------------------
+  // This function unloads MFT clusters
+  //--------------------------------------------------------------------
+
+  for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
+    fMFTClusterArray[iPlane]      -> Clear("C");
+    fMFTClusterArrayFront[iPlane] -> Clear("C");
+    fMFTClusterArrayBack[iPlane]  -> Clear("C");
+  }
+
+}
+
+//====================================================================================================================================================
+
+Int_t AliMFTTrackerMU::Clusters2Tracks(AliESDEvent *event) {
+
+  //--------------------------------------------------------------------
+  // This functions reconstructs the Muon Forward Tracks
+  // The clusters must be already loaded !
+  //--------------------------------------------------------------------
+
+  TTree *outputTreeMuonGlobalTracks = new TTree("AliMuonForwardTracks", "Tree of AliMuonForwardTracks");
+  TClonesArray *muonForwardTracks = new TClonesArray("AliMuonForwardTrack");
+  outputTreeMuonGlobalTracks -> Branch("tracks", &muonForwardTracks);
+  //--------------------------------------------------------------------
+
+  fESD = event;
+
+  //----------- Read ESD MUON tracks -------------------
+
+  Int_t nTracksMUON = event->GetNumberOfMuonTracks();
+
+  AliInfo(Form("Number of ESD MUON tracks: %d\n", nTracksMUON));
+
+  Int_t iTrack=0;
+  while (iTrack<nTracksMUON) {
+
+    fNPlanesMFTAnalyzed = 0;
+
+    AliDebug(1, "**************************************************************************************\n");
+    AliDebug(1, Form("***************************   MUON TRACK %3d   ***************************************\n", iTrack));
+    AliDebug(1, "**************************************************************************************\n");
+    
+    fCandidateTracks -> Delete();
+    
+    fNPlanesMFTAnalyzed = 0;
+    
+    const AliESDMuonTrack *esdTrack = event->GetMuonTrack(iTrack);
+    fMUONTrack = NULL;
+    AliMUONESDInterface::ESDToMUON(*esdTrack, *fMUONTrack, kFALSE);
+
+    // the track we are going to build, starting from fMUONTrack and adding the MFT clusters
+    AliMuonForwardTrack *track = new ((*fCandidateTracks)[0]) AliMuonForwardTrack();
+    track -> SetMUONTrack(new AliMUONTrack(*fMUONTrack));
+    track -> SetMCLabel(fMUONTrack->GetMCLabel());
+    track -> SetMatchTrigger(fMUONTrack->GetMatchTrigger());
+
+    //------------------------- NOW THE CYCLE OVER THE MFT PLANES STARTS ---------------------------------------
+
+    for (Int_t iPlane=fNPlanesMFT-1; iPlane>=0; iPlane--) {   /* *** do not reverse the order of this cycle!!! 
+                                                                *** this reflects the fact that the extrapolation is performed 
+                                                                *** starting from the last MFT plane back to the origin */
+      
+      // --------- updating the array of candidates according to the clusters available in the i-th plane ---------
+      
+      fNPlanesMFTAnalyzed++;
+      
+      Int_t nCandidates = fCandidateTracks->GetEntriesFast();
+      for (Int_t iCandidate=0; iCandidate<nCandidates; iCandidate++) {
+       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 
+       // (several new tracks can be created for one old track)
+       if (FindClusterInPlane(iPlane) == kDiverged) {
+         fGlobalTrackingDiverged = kTRUE;
+         break;
+       }
+
+       if ((fNPlanesMFTAnalyzed-fCurrentTrack->GetNMFTClusters())>fNMaxMissingMFTClusters || fIsPlaneMandatory[iPlane]) {
+         fCandidateTracks->Remove(fCurrentTrack);     // the old track is removed after the check;
+       }
+      }
+      if (fGlobalTrackingDiverged) {
+       if (fScaleSigmaClusterCut>0) fScaleSigmaClusterCut -= 0.1;
+       continue;
+      }
+
+      fCandidateTracks->Compress();
+      
+    }      
+    
+    // -------------------------- END OF THE CYCLE OVER THE MFT PLANES --------------------------------------------
+    
+    fGlobalTrackingDiverged = kFALSE;
+    fScaleSigmaClusterCut = 1.0;
+    
+    AliDebug(1, "Finished cycle over planes");
+    
+    iTrack++;
+
+    // If we have several final tracks, we must find the best candidate:
+    
+    Int_t nFinalTracks = fCandidateTracks->GetEntriesFast();
+    AliDebug(1, Form("nFinalTracks = %d", nFinalTracks));
+    
+    Int_t nGoodClustersBestCandidate =  0;
+    Int_t idBestCandidate            =  0;
+    Double_t bestChi2                = -1.;  // variable defining the best candidate
+
+    for (Int_t iFinalCandidate=0; iFinalCandidate<nFinalTracks; iFinalCandidate++) {
+      
+      AliMuonForwardTrack *finalTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iFinalCandidate);
+      Int_t nMFTClusters  = finalTrack->GetNMFTClusters();
+
+      Double_t chi2 = 0;
+      for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) {
+       AliMFTCluster *localCluster = finalTrack->GetMFTCluster(iCluster);
+        chi2 += localCluster->GetLocalChi2();
+      }
+      chi2 /= nMFTClusters;
+
+      // now comparing the tracks in order to find the best one
+      
+      if (chi2<bestChi2 || bestChi2<0) {
+       bestChi2 = chi2;
+       idBestCandidate = iFinalCandidate;
+      }
+      
+    }
+    
+    if (nFinalTracks) {
+
+      AliMuonForwardTrack *newTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(idBestCandidate);
+      newTrack -> SetNWrongClustersMC(newTrack->GetNMFTClusters() - nGoodClustersBestCandidate);
+
+      //----------------------- Save the information to the AliESDMuonForwardTrack object
+
+      AliESDMuonGlobalTrack *myESDTrack = event->NewMuonGlobalTrack();
+      myESDTrack -> SetPxPyPz(newTrack->Px(), newTrack->Py(), newTrack->Pz());
+      myESDTrack -> SetChi2(newTrack->GetGlobalChi2());
+      myESDTrack -> SetCharge(newTrack->GetCharge());
+      myESDTrack -> SetMatchTrigger(newTrack->GetMatchTrigger());
+      
+      //---------------------------------------------------------------------------------
+
+      new ((*muonForwardTracks)[muonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack);
+
+    }
+
+    fCandidateTracks->Delete();
+    fFinalBestCandidate = NULL;
+   
+  }
+
+  Int_t myEventID = 0;
+  TFile *outputFileMuonGlobalTracks = new TFile("MuonGlobalTracks.root", "update");
+  while (outputFileMuonGlobalTracks->cd(Form("Event%d",myEventID))) myEventID++;
+  outputFileMuonGlobalTracks -> mkdir(Form("Event%d",myEventID));
+  outputFileMuonGlobalTracks -> cd(Form("Event%d",myEventID));
+  outputTreeMuonGlobalTracks -> Write();
+  outputFileMuonGlobalTracks -> Close();
+
+  return 0;
+
+}
+
+//=========================================================================================================================================
+
+void AliMFTTrackerMU::SeparateFrontBackClusters() {
+
+  for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
+    fMFTClusterArrayFront[iPlane]->Delete();
+    fMFTClusterArrayBack[iPlane] ->Delete();
+    for (Int_t iCluster=0; iCluster<fMFTClusterArray[iPlane]->GetEntries(); iCluster++) {
+      AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->At(iCluster);
+      if (TMath::Abs(cluster->GetZ())<TMath::Abs(fSegmentation->GetPlane(iPlane)->GetZCenter())) {
+       new ((*fMFTClusterArrayFront[iPlane])[fMFTClusterArrayFront[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
+      }
+      else {
+       new ((*fMFTClusterArrayBack[iPlane])[fMFTClusterArrayBack[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
+      }
+    }
+  }
+
+}
+
+//==========================================================================================================================================
+
+Int_t AliMFTTrackerMU::FindClusterInPlane(Int_t planeId) { 
+  
+  AliDebug(2, Form(">>>> executing AliMuonForwardTrackFinder::FindClusterInPlane(%d)\n", planeId));
+
+  // !!!!!!!!! coordinates and errors on the interaction vertex should be taken from the event itself (ITS) if available
+
+  // propagate track to plane #planeId (both to front and back active sensors)
+  // look for compatible clusters
+  // update TrackParam at found cluster (if any) using Kalman Filter
+
+  AliMUONTrackParam currentParamFront, currentParamBack, currentParamForResearchFront, currentParamForResearchBack;
+
+  if (planeId == fNPlanesMFT-1) {      // last plane of the telecope
+    currentParamFront = (*((AliMUONTrackParam*)(fMUONTrack->GetTrackParamAtCluster()->First())));
+    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); 
+    }
+    else {
+      AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamFront, zExtrap);
+      AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamBack,  zExtrap);
+    }
+    AliMUONTrackExtrap::ExtrapToVertex(&currentParamForResearchFront, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY); 
+    AliMUONTrackExtrap::ExtrapToVertex(&currentParamForResearchBack,  xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY); 
+  }
+  else {          // MFT planes others than the last one: mult. scattering correction because of the upstream MFT planes is performed
+    currentParamFront = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
+    currentParamBack  = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
+    currentParamForResearchFront = currentParamFront;
+    currentParamForResearchBack  = currentParamBack;
+    AliMUONTrackExtrap::AddMCSEffect(&currentParamFront,           (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
+                                                                   fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
+    AliMUONTrackExtrap::AddMCSEffect(&currentParamForResearchFront,(fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
+                                                                   fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
+    AliMUONTrackExtrap::AddMCSEffect(&currentParamBack,            (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
+                                                                   fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
+    AliMUONTrackExtrap::AddMCSEffect(&currentParamForResearchBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
+                                                                   fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
+  }
+  // for all planes: extrapolation to the Z of the plane
+  AliMUONTrackExtrap::ExtrapToZCov(&currentParamFront,            -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());   
+  AliMUONTrackExtrap::ExtrapToZCov(&currentParamForResearchFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());
+  AliMUONTrackExtrap::ExtrapToZCov(&currentParamBack,             -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());   
+  AliMUONTrackExtrap::ExtrapToZCov(&currentParamForResearchBack,  -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());
+
+  //---------------------------------------------------------------------------------------
+
+  TMatrixD covFront(5,5); covFront = currentParamForResearchFront.GetCovariances();
+  TMatrixD covBack(5,5);  covBack  = currentParamForResearchBack.GetCovariances();
+  
+  Double_t squaredError_X_Front = covFront(0,0);
+  Double_t squaredError_Y_Front = covFront(2,2);
+  Double_t squaredError_X_Back  = covBack(0,0);
+  Double_t squaredError_Y_Back  = covBack(2,2);
+
+  Double_t corrFact = 1.0;
+
+  Double_t researchRadiusFront = TMath::Sqrt(squaredError_X_Front + squaredError_Y_Front);
+  Double_t researchRadiusBack  = TMath::Sqrt(squaredError_X_Back  + squaredError_Y_Back);
+  if (0.5*(researchRadiusFront+researchRadiusBack)<fMinResearchRadiusAtPlane[planeId]) {
+    corrFact = fMinResearchRadiusAtPlane[planeId]/(0.5*(researchRadiusFront+researchRadiusBack));
+  }
+
+  //---------------------------------------------------------------------------------------
+
+  Double_t chi2cut = 2.*fScaleSigmaClusterCut*fScaleSigmaClusterCut*fSigmaClusterCut*fSigmaClusterCut;     // depends on the number of variables (here, 2)
+  
+  // Analyizing the clusters: FRONT ACTIVE ELEMENTS
+  
+  Int_t nClustersFront = fMFTClusterArrayFront[planeId]->GetEntries();
+  AliDebug(2, Form("There are %3d clusters in plane %02d FRONT\n", nClustersFront, planeId));
+  
+  for (Int_t iCluster=0; iCluster<nClustersFront; iCluster++) {
+
+    Bool_t isGoodChi2 = kFALSE;
+
+    AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayFront[planeId]->At(iCluster); 
+    Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchFront, cluster);     // describes the compatibility between the track and the cluster
+    if (chi2<chi2cut) isGoodChi2 = kTRUE;
+
+    if (isGoodChi2) {
+      AliDebug(3, Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
+      AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
+      if (fCandidateTracks->GetEntriesFast() > fMaxNCandidates) return kDiverged;
+      newTrack->AddTrackParamAtMFTCluster(currentParamFront, *cluster);    // creating new track param and attaching the cluster
+      AliDebug(2, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)", 
+                      planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters()));
+      newTrack->SetPlaneExists(planeId);
+    }
+    else AliDebug(3, Form("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
+
+  }
+
+  // Analyizing the clusters: BACK ACTIVE ELEMENTS
+  
+  Int_t nClustersBack = fMFTClusterArrayBack[planeId]->GetEntries();
+  AliDebug(2, Form("There are %3d clusters in plane %02d BACK\n", nClustersBack, planeId));
+  
+  for (Int_t iCluster=0; iCluster<nClustersBack; iCluster++) {
+
+    Bool_t isGoodChi2 = kFALSE;
+
+    AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayBack[planeId]->At(iCluster); 
+    Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchBack, cluster);     // describes the compatibility between the track and the cluster
+    if (chi2<chi2cut) isGoodChi2 = kTRUE;
+
+    if (isGoodChi2) {
+      AliDebug(3,Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
+      AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
+      if (fCandidateTracks->GetEntriesFast() > fMaxNCandidates) return kDiverged;
+      newTrack->AddTrackParamAtMFTCluster(currentParamBack, *cluster);    // creating new track param and attaching the cluster
+      AliDebug(2, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)", 
+                      planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters()));
+      newTrack->SetPlaneExists(planeId);
+    }
+    else AliDebug(3,Form("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
+
+  }
+
+  //---------------------------------------------------------------------------------------------
+
+  return kConverged;
+  
+}
+
+//==========================================================================================================================================
+
+Double_t AliMFTTrackerMU::TryOneCluster(const AliMUONTrackParam &trackParam, AliMFTCluster *cluster) {
+
+  // Test the compatibility between the track and the cluster (using trackParam's covariance matrix):
+  // return the corresponding Chi2
+  // assume the track parameters are given at the Z of the cluster
+  
+  // Set differences between trackParam and cluster in the bending and non bending directions
+  Double_t dX = cluster->GetX() - trackParam.GetNonBendingCoor();
+  Double_t dY = cluster->GetY() - trackParam.GetBendingCoor();
+  AliDebug(3,Form("dX = %f, dY = %f\n", dX, dY));
+  
+  // Calculate errors and covariances
+  const TMatrixD& kParamCov = trackParam.GetCovariances();
+  Double_t sigmaX2 = kParamCov(0,0) + cluster->GetErrX2();
+  Double_t sigmaY2 = kParamCov(2,2) + cluster->GetErrY2();
+  AliDebug(3, Form("dX2 = %f, dY2 = %f\n", sigmaX2, sigmaY2));
+  Double_t covXY   = kParamCov(0,2);
+  Double_t det     = sigmaX2 * sigmaY2 - covXY * covXY;
+  
+  // Compute chi2
+  if (det==0.) return 1.e10;
+  return (dX*dX*sigmaY2 + dY*dY*sigmaX2 - 2.*dX*dY*covXY) / det;
+  
+}
+
+//=========================================================================================================================================
+
diff --git a/MFT/AliMFTTrackerMU.h b/MFT/AliMFTTrackerMU.h
new file mode 100644 (file)
index 0000000..b1703f9
--- /dev/null
@@ -0,0 +1,116 @@
+#ifndef AliMFTTrackerMU_H
+#define AliMFTTrackerMU_H
+
+/* Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//====================================================================================================================================================
+//
+//                          MFT tracker
+//
+// Class for the creation of the "global muon tracks" built from the tracks reconstructed in the 
+// muon spectrometer and the clusters of the Muon Forward Tracker
+//
+// Contact author: antonio.uras@cern.ch
+//
+//====================================================================================================================================================
+
+#include "TTree.h"
+#include "AliLog.h"
+#include "AliGeomManager.h"
+#include "AliESDEvent.h"
+#include "AliESDMuonTrack.h"
+#include "AliESDMuonGlobalTrack.h"
+#include "AliTracker.h"
+#include "AliRun.h"
+#include "AliMFT.h"
+#include "AliMUONTrackExtrap.h"
+#include "AliMUONTrack.h"
+#include "AliMUONESDInterface.h"
+#include "AliMuonForwardTrack.h"
+
+//====================================================================================================================================================
+
+class AliMFTTrackerMU : public AliTracker {
+
+public:
+
+  enum {kConverged, kDiverged};
+
+  AliMFTTrackerMU();
+  virtual ~AliMFTTrackerMU();
+
+  Int_t LoadClusters(TTree *cf);
+  void UnloadClusters();
+  Int_t Clusters2Tracks(AliESDEvent *event);
+
+  void SetNPlanesMFT(Int_t nPlanesMFT) { fNPlanesMFT = nPlanesMFT; }
+  void SeparateFrontBackClusters();
+
+  Int_t FindClusterInPlane(Int_t planeId);
+
+  void SetVertexError(Double_t xErr, Double_t yErr, Double_t zErr) { fVertexErrorX=xErr; fVertexErrorY=yErr; fVertexErrorZ=zErr; }
+  void SetMinResearchRadiusAtPlane(Int_t plane, Double_t radius) { if (plane>=0 && plane<fNMaxPlanes) fMinResearchRadiusAtPlane[plane] = radius; }
+
+  Double_t TryOneCluster(const AliMUONTrackParam &trackParam, AliMFTCluster *cluster);
+
+  /// Dummy implementation
+  virtual Int_t PropagateBack(AliESDEvent* /*event*/) {return 0;}
+  /// Dummy implementation
+  virtual Int_t RefitInward(AliESDEvent* /*event*/) {return 0;}
+  /// Dummy implementation
+  virtual AliCluster *GetCluster(Int_t /*index*/) const {return 0;}
+
+
+protected:
+
+  static const Int_t fNMaxPlanes = AliMFTConstants::fNMaxPlanes;        // max number of MFT planes
+  static const Double_t fRadLengthSi;
+  static const Int_t fMaxNCandidates = 1000;
+
+  AliESDEvent *fESD;   //! pointer to the ESD event
+
+  AliMFT *fMFT;                        //!
+  AliMFTSegmentation *fSegmentation;   //!
+  Int_t fNPlanesMFT, fNPlanesMFTAnalyzed;
+
+  Double_t fSigmaClusterCut;         // to select the clusters in the MFT planes which are compatible with the extrapolated muon track
+  Double_t fScaleSigmaClusterCut;    // to tune the cut on the compatible clusters in case of too many candidates
+
+  Int_t fNMaxMissingMFTClusters;              // max. number of MFT clusters which can be missed in the global fit procedure
+  Bool_t fIsPlaneMandatory[fNMaxPlanes];      // specifies which MFT planes cannot be missed in the global fit procedure
+
+  Bool_t fGlobalTrackingDiverged;    // to keep memory of a possible divergence in the global tracking finding
+
+  TClonesArray *fMFTClusterArray[fNMaxPlanes];         //! array of clusters for the planes of the MFT
+  TClonesArray *fMFTClusterArrayFront[fNMaxPlanes];    //! array of front clusters for the planes of the MFT
+  TClonesArray *fMFTClusterArrayBack[fNMaxPlanes];     //! array of back clusters for the planes of the MFT
+
+  TClonesArray *fCandidateTracks;   //! array of candidate global tracks 
+
+  AliMUONTrack *fMUONTrack;  //! muon track being analyzed
+
+  AliMuonForwardTrack *fCurrentTrack;     //! muon extrapolated track being tested
+  AliMuonForwardTrack *fFinalBestCandidate;     //! best final candidate (if any)
+
+  // uncertainty on the x position of the primary vertex (seed for the extrapolation of the MUON tracks)
+  Double_t fVertexErrorX;   
+  Double_t fVertexErrorY;
+  Double_t fVertexErrorZ;
+
+  Bool_t fBransonCorrection;    // if TRUE, Branson Correction is applied when extrapolating the MUON tracks to the vertex region
+
+  Double_t fMinResearchRadiusAtPlane[fNMaxPlanes];
+
+private:
+
+  AliMFTTrackerMU(const AliMFTTrackerMU &tracker);
+  AliMFTTrackerMU & operator=(const AliMFTTrackerMU &tracker);
+  ClassDef(AliMFTTrackerMU,1)   //MFT tracker for muon tracks
+
+};
+
+//====================================================================================================================================================
+
+#endif 
+
index 4d39448..9ed3236 100644 (file)
@@ -576,7 +576,7 @@ Int_t AliMuonForwardTrackFinder::LoadNextTrack() {
        fCurrentTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iTrack);
        // if the old track is compatible with the new cluster, the track is updated and inserted as new track in the array 
        // (several new tracks can be created for one old track)
-       if (FindClusterInPlane(iPlane) == 2) {
+       if (FindClusterInPlane(iPlane) == kDiverged) {
          fGlobalTrackingDiverged = kTRUE;
          if (fScaleSigmaClusterCut>0) fScaleSigmaClusterCut -= 0.1;
          return 6;
@@ -927,7 +927,7 @@ Int_t AliMuonForwardTrackFinder::FindClusterInPlane(Int_t planeId) {
     if (isGoodChi2) {
       AliDebug(3, Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
       AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
-      if (fCandidateTracks->GetEntriesFast() > fMaxNCandidates) return 2;
+      if (fCandidateTracks->GetEntriesFast() > fMaxNCandidates) return kDiverged;
       newTrack->AddTrackParamAtMFTCluster(currentParamFront, *cluster);    // creating new track param and attaching the cluster
       AliDebug(2, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)", 
                       planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters()));
@@ -981,7 +981,7 @@ Int_t AliMuonForwardTrackFinder::FindClusterInPlane(Int_t planeId) {
     if (isGoodChi2) {
       AliDebug(3,Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
       AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
-      if (fCandidateTracks->GetEntriesFast() > fMaxNCandidates) return 2;
+      if (fCandidateTracks->GetEntriesFast() > fMaxNCandidates) return kDiverged;
       newTrack->AddTrackParamAtMFTCluster(currentParamBack, *cluster);    // creating new track param and attaching the cluster
       AliDebug(2, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)", 
                       planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters()));
@@ -1013,7 +1013,7 @@ Int_t AliMuonForwardTrackFinder::FindClusterInPlane(Int_t planeId) {
     }
   }
 
-  return 0;
+  return kConverged;
   
 }
 
index b69dfa7..6f0e17d 100644 (file)
@@ -63,6 +63,8 @@ class AliMuonForwardTrackFinder : public TObject {
   
 public:
 
+  enum {kConverged, kDiverged};
+
   enum matchingOption {kRealMatching, kIdealMatching};
 
   AliMuonForwardTrackFinder();
index 58c302c..64bf1ae 100644 (file)
@@ -25,7 +25,7 @@
 # SHLIBS - Shared Libraries and objects for linking (Executables only)           #
 #--------------------------------------------------------------------------------#
 
-set ( SRCS AliMFTReconstructor.cxx AliMFTClusterFinder.cxx AliMuonForwardTrack.cxx AliMuonForwardTrackPair.cxx AliESDEventMFT.cxx AliMuonForwardTrackFinder.cxx AliMFTClusterQA.cxx AliMFTRecoParam.cxx )
+set ( SRCS AliMFTReconstructor.cxx AliMFTTrackerMU.cxx AliMFTClusterFinder.cxx AliMuonForwardTrack.cxx AliMuonForwardTrackPair.cxx AliESDEventMFT.cxx AliMuonForwardTrackFinder.cxx AliMFTClusterQA.cxx AliMFTRecoParam.cxx )
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
 
index d939e3f..a3cbe39 100644 (file)
@@ -7,6 +7,7 @@
 #pragma link off all functions;
  
 #pragma link C++ class  AliMFTReconstructor+;
+#pragma link C++ class  AliMFTTrackerMU+;
 #pragma link C++ class  AliMFTClusterFinder+;
 #pragma link C++ class  AliMuonForwardTrack+;
 #pragma link C++ class  AliMuonForwardTrackPair+;