Updated MFT version, some problems from the previous attempt are fixed (Antonio)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 2 Feb 2012 12:06:32 +0000 (12:06 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 2 Feb 2012 12:06:32 +0000 (12:06 +0000)
31 files changed:
MFT/AliMFT.cxx
MFT/AliMFT.h
MFT/AliMFTCluster.cxx
MFT/AliMFTCluster.h
MFT/AliMFTClusterFinder.cxx
MFT/AliMFTClusterFinder.h
MFT/AliMFTClusterQA.cxx
MFT/AliMFTClusterQA.h
MFT/AliMFTConstants.cxx [new file with mode: 0644]
MFT/AliMFTConstants.h [new file with mode: 0644]
MFT/AliMFTDigit.cxx
MFT/AliMFTDigit.h
MFT/AliMFTDigitizer.cxx
MFT/AliMFTDigitizer.h
MFT/AliMFTGeometry.root
MFT/AliMFTPlane.cxx
MFT/AliMFTPlane.h
MFT/AliMFTReconstructor.cxx
MFT/AliMFTReconstructor.h
MFT/AliMFTSegmentation.cxx
MFT/AliMFTSegmentation.h
MFT/AliMuonForwardTrack.cxx
MFT/AliMuonForwardTrack.h
MFT/AliMuonForwardTrackFinder.C
MFT/AliMuonForwardTrackFinder.cxx
MFT/AliMuonForwardTrackFinder.h
MFT/AliMuonForwardTrackPair.cxx
MFT/AliMuonForwardTrackPair.h
MFT/CMakelibMFTbase.pkg
MFT/MFTbaseLinkDef.h
MFT/SetMFTGeometry.C

index 1a23486..f4b9a4c 100644 (file)
@@ -63,7 +63,7 @@ AliMFT::AliMFT():
   fChargeDispersion(0),
   fSingleStepForChargeDispersion(0),
   fNStepForChargeDispersion(0),
-  fDensitySiOverSupport(10)
+  fDensitySupportOverSi(0.1)
 {
 
   // default constructor
@@ -86,7 +86,7 @@ AliMFT::AliMFT(const Char_t *name, const Char_t *title):
   fChargeDispersion(0),
   fSingleStepForChargeDispersion(0),
   fNStepForChargeDispersion(0),
-  fDensitySiOverSupport(10)
+  fDensitySupportOverSi(0.1)
 {
 
   fNameGeomFile = "AliMFTGeometry.root";
@@ -113,7 +113,7 @@ AliMFT::AliMFT(const Char_t *name, const Char_t *title, Char_t *nameGeomFile):
   fChargeDispersion(0),
   fSingleStepForChargeDispersion(0),
   fNStepForChargeDispersion(0),
-  fDensitySiOverSupport(10)
+  fDensitySupportOverSi(0.1)
 {
 
   fNameGeomFile = nameGeomFile;
@@ -175,7 +175,7 @@ void AliMFT::CreateMaterials() {
   AliMaterial(++matId, "Readout", aSi,   zSi,    dSi,    radSi,  absSi  );  
   AliMedium(kReadout,  "Readout", matId, unsens, isxfld, sxmgmx, tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
 
-  AliMaterial(++matId, "Support", aSi,   zSi,    dSi/fDensitySiOverSupport, fDensitySiOverSupport*radSi, fDensitySiOverSupport*absSi);  
+  AliMaterial(++matId, "Support", aSi,   zSi,    dSi*fDensitySupportOverSi, radSi/fDensitySupportOverSi, absSi/fDensitySupportOverSi);  
   AliMedium(kSupport,  "Support", matId, unsens, isxfld,  sxmgmx,    tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
     
   AliInfo("End MFT materials");
@@ -205,6 +205,8 @@ void AliMFT::StepManager() {
 
   // Full Step Manager
 
+  AliDebug(2, Form("Entering StepManager: gMC->CurrentVolName() = %s", gMC->CurrentVolName()));
+
   if (!fSegmentation) AliFatal("No segmentation available");    // DO WE HAVE A SEGMENTATION???
 
   if (!(this->IsActive())) return;
@@ -247,8 +249,8 @@ void AliMFT::StepManager() {
   gMC->TrackPosition(position);
   gMC->TrackMomentum(momentum);
 
-  AliDebug(1, Form("AliMFT::StepManager()->%s Hit #%06d (z=%f) belongs to track %02d\n", 
-                  gMC->CurrentVolName(), fNhits, position.Z(), gAlice->GetMCApp()->GetCurrentTrackNumber())); 
+  AliDebug(1, Form("AliMFT::StepManager()->%s Hit #%06d (x=%f, y=%f, z=%f) belongs to track %02d\n", 
+                  gMC->CurrentVolName(), fNhits, position.X(), position.Y(), position.Z(), gAlice->GetMCApp()->GetCurrentTrackNumber())); 
 
   hit.SetPosition(position);
   hit.SetTOF(gMC->TrackTime());
@@ -309,6 +311,8 @@ TGeoVolumeAssembly* AliMFT::CreateVol() {
                                                    plane->GetRMaxSupport(),
                                                    0.5*(plane->GetSupportElement(0)->GetAxis(2)->GetXmax() - 
                                                         plane->GetSupportElement(0)->GetAxis(2)->GetXmin()) );
+    AliDebug(2, Form("Created vol %s", supportElem->GetName()));
+    supportElem->SetLineColor(kCyan-9);
     vol -> AddNode(supportElem, 0, new TGeoTranslation(origin[0], origin[1], origin[2]));
 
     AliDebug(1, "support elements created!");
@@ -328,6 +332,8 @@ TGeoVolumeAssembly* AliMFT::CreateVol() {
       for (Int_t iSlice=0; iSlice<fNSlices; iSlice++) {
        origin[2] = -0.5*(plane->GetActiveElement(iActive)->GetAxis(2)->GetXmin() + 2*dz*(iSlice+1) + plane->GetActiveElement(iActive)->GetAxis(2)->GetXmin() + 2*dz*(iSlice) );
        TGeoVolume *activeElem = gGeoManager->MakeBox(Form("MFT_plane%02d_active%03d_slice%02d", iPlane, iActive, iSlice), silicon, dx, dy, dz);
+       AliDebug(2, Form("Created vol %s", activeElem->GetName()));
+       activeElem->SetLineColor(kGreen);
        vol -> AddNode(activeElem, 0, new TGeoTranslation(origin[0], origin[1], origin[2]));
       }
 
@@ -348,6 +354,8 @@ TGeoVolumeAssembly* AliMFT::CreateVol() {
       origin[2] = -0.5*(plane->GetReadoutElement(iReadout)->GetAxis(2)->GetXmax() + plane->GetReadoutElement(iReadout)->GetAxis(2)->GetXmin());
 
       TGeoVolume *readoutElem = gGeoManager->MakeBox(Form("MFT_plane%02d_readout%03d", iPlane, iReadout), readout, dx, dy, dz);
+      AliDebug(2, Form("Created vol %s", readoutElem->GetName()));
+      readoutElem->SetLineColor(kRed);
       vol -> AddNode(readoutElem, 0, new TGeoTranslation(origin[0], origin[1], origin[2]));
       
     }
index c8d1c3e..884574c 100644 (file)
@@ -24,7 +24,6 @@
 #include "AliDetector.h"
 #include "AliMC.h"
 #include "AliMagF.h"
-#include "AliMFT.h"
 #include "AliMFTHit.h"
 #include "AliMFTDigit.h"
 #include "AliMFTCluster.h"
@@ -34,6 +33,7 @@
 #include "AliMFTPlane.h"
 #include "TString.h"
 #include "TObjArray.h"
+#include "AliMFTConstants.h"
 
 //====================================================================================================================================================
 
@@ -64,20 +64,20 @@ public:
   void CreateDigits();
   void CreateRecPoints();
 
-  TObjArray*    GetSDigitsList()            const { return fSDigitsPerPlane; }                                                  // get sdigits list for all planes
-  TClonesArray* GetSDigitsList(Int_t plane) const { return fSDigitsPerPlane ? (TClonesArray*) fSDigitsPerPlane->At(plane):0; }  // get sdigits list for a plane
+  TObjArray*    GetSDigitsList()            const { return fSDigitsPerPlane; }     // get sdigits list for all planes
+  TClonesArray* GetSDigitsList(Int_t plane) const { return fSDigitsPerPlane ? (TClonesArray*) fSDigitsPerPlane->At(plane):0; } 
 
-  TObjArray*    GetDigitsList()            const{return fDigitsPerPlane;}                                                   // get digits list for all layers
-  TClonesArray* GetDigitsList(Int_t plane) const{return fDigitsPerPlane ? (TClonesArray*) fDigitsPerPlane->At(plane):0; }   // get digits list for a plane
+  TObjArray*    GetDigitsList()            const{return fDigitsPerPlane;}          // get digits list for all layers
+  TClonesArray* GetDigitsList(Int_t plane) const{return fDigitsPerPlane ? (TClonesArray*) fDigitsPerPlane->At(plane):0; }
 
-  TObjArray*    GetRecPointsList()            const{return fRecPointsPerPlane;}                                                      // get digits list for all layers
-  TClonesArray* GetRecPointsList(Int_t plane) const{return fRecPointsPerPlane ? (TClonesArray*) fRecPointsPerPlane->At(plane):0; }   // get digits list for a plane
+  TObjArray*    GetRecPointsList()            const{return fRecPointsPerPlane;}    // get cluster list for all layers
+  TClonesArray* GetRecPointsList(Int_t plane) const{return fRecPointsPerPlane ? (TClonesArray*) fRecPointsPerPlane->At(plane):0; }
 
   void ResetSDigits()   { if(fSDigitsPerPlane)   for(int iPlane=0; iPlane<fNPlanes; iPlane++) ((TClonesArray*) fSDigitsPerPlane  ->At(iPlane))->Clear(); }   // reset sdigits list  
   void ResetDigits()    { if(fDigitsPerPlane)    for(int iPlane=0; iPlane<fNPlanes; iPlane++) ((TClonesArray*) fDigitsPerPlane   ->At(iPlane))->Clear(); }   // reset digits list
   void ResetRecPoints() { if(fRecPointsPerPlane) for(int iPlane=0; iPlane<fNPlanes; iPlane++) ((TClonesArray*) fRecPointsPerPlane->At(iPlane))->Clear(); }   // reset recPoints list
   
-  AliDigitizer* CreateDigitizer(AliDigitizationInput *digInp) const { return new AliMFTDigitizer(digInp); }   // from AliModule invoked from AliSimulation::RunDigitization()
+  AliDigitizer* CreateDigitizer(AliDigitizationInput *digInp) const { return new AliMFTDigitizer(digInp); }
   
   AliMFTSegmentation* GetSegmentation() const { return fSegmentation; }
 
@@ -99,11 +99,11 @@ public:
   Int_t GetNStepForChargeDispersion() { return fNStepForChargeDispersion; }
   Double_t GetSingleStepForChargeDispersion() { return fSingleStepForChargeDispersion; }
 
-  void SetDensitySiOverSupport(Double_t density) { fDensitySiOverSupport = density; }
+  void SetDensitySupportOverSi(Double_t density) { if (density>1e-6) fDensitySupportOverSi=density; else fDensitySupportOverSi=1e-6; }
 
 protected:
 
-  static const Int_t fNMaxPlanes = 20;        // max number of MFT planes
+  static const Int_t fNMaxPlanes = AliMFTConstants::fNMaxPlanes;        // max number of MFT planes
 
   Int_t fVersion;
 
@@ -123,7 +123,7 @@ protected:
   Double_t fSingleStepForChargeDispersion;
   Int_t fNStepForChargeDispersion;
 
-  Double_t fDensitySiOverSupport;
+  Double_t fDensitySupportOverSi;
   
 private:
 
index 3cfb7a2..1870894 100644 (file)
@@ -24,6 +24,8 @@
 #include "TObject.h"
 #include "AliMUONRawCluster.h"
 #include "AliMUONVCluster.h"
+#include "AliMFTDigit.h"
+#include "TMath.h"
 #include "AliMFTCluster.h"
 
 ClassImp(AliMFTCluster)
@@ -40,16 +42,175 @@ AliMFTCluster::AliMFTCluster():
   fErrZ(0),
   fNElectrons(0),
   fNMCTracks(0),
-  fPlane(0),
+  fPlane(-1),
+  fDetElemID(-1),
   fSize(0),
   fTrackChi2(0),
-  fLocalChi2(0)
+  fLocalChi2(0),
+  fDigitsInCluster(0),
+  fIsClusterEditable(kTRUE)
 {
 
   // default constructor
 
   for (Int_t iTrack=0; iTrack<fNMaxMCTracks; iTrack++) fMCLabel[iTrack] = -1;
 
+  fDigitsInCluster = new TClonesArray("AliMFTDigit", fNMaxDigitsPerCluster);
+
+}
+
+//====================================================================================================================================================
+
+AliMFTCluster::AliMFTCluster(const AliMFTCluster& cluster): 
+  TObject(cluster),
+  fX(cluster.fX), 
+  fY(cluster.fY), 
+  fZ(cluster.fZ),
+  fErrX(cluster.fErrX), 
+  fErrY(cluster.fErrY), 
+  fErrZ(cluster.fErrZ),
+  fNElectrons(cluster.fNElectrons),
+  fNMCTracks(cluster.fNMCTracks),
+  fPlane(cluster.fPlane),
+  fDetElemID(cluster.fDetElemID),
+  fSize(cluster.fSize),
+  fTrackChi2(cluster.fTrackChi2),
+  fLocalChi2(cluster.fLocalChi2),
+  fDigitsInCluster(cluster.fDigitsInCluster),
+  fIsClusterEditable(cluster.fIsClusterEditable)
+{
+
+  // copy constructor
+  for (Int_t iTrack=0; iTrack<fNMaxMCTracks; iTrack++) fMCLabel[iTrack] = (cluster.fMCLabel)[iTrack];
+  
+}
+
+//====================================================================================================================================================
+
+AliMFTCluster& AliMFTCluster::operator=(const AliMFTCluster& cluster) {
+
+  // Asignment operator
+
+  // check assignement to self
+  if (this == &cluster) return *this;
+
+  // base class assignement
+  TObject::operator=(cluster);
+  
+  // clear memory
+  Clear();
+  
+  fX                 = cluster.fX; 
+  fY                 = cluster.fY; 
+  fZ                 = cluster.fZ;
+  fErrX              = cluster.fErrX; 
+  fErrY              = cluster.fErrY; 
+  fErrZ              = cluster.fErrZ;
+  fNElectrons        = cluster.fNElectrons;
+  fNMCTracks         = cluster.fNMCTracks;
+  fPlane             = cluster.fPlane;
+  fDetElemID         = cluster.fDetElemID;
+  fSize              = cluster.fSize;
+  fTrackChi2         = cluster.fTrackChi2;
+  fLocalChi2         = cluster.fLocalChi2;
+  fDigitsInCluster   = cluster.fDigitsInCluster;
+  fIsClusterEditable = cluster.fIsClusterEditable;
+
+  for (Int_t iTrack=0; iTrack<fNMaxMCTracks; iTrack++) fMCLabel[iTrack] = (cluster.fMCLabel)[iTrack];
+  
+  return *this;
+
+}
+
+//====================================================================================================================================================
+
+Double_t AliMFTCluster::GetDistanceFromPixel(AliMFTDigit *pixel) {
+
+  // the distance is expressed in units of pixels!!!
+  // useful to decide if the pixel is compatible with the current digits array
+
+  Double_t distance = -1;
+
+  if (!fSize) return distance;
+
+  if (pixel->GetDetElemID()!=fDetElemID || pixel->GetPlane()!=fPlane) return 9999.;
+
+  for (Int_t iDigit=0; iDigit<fDigitsInCluster->GetEntries(); iDigit++) {
+    AliMFTDigit *tmpDig = (AliMFTDigit*) fDigitsInCluster->At(iDigit);
+    Int_t distX = TMath::Abs(tmpDig->GetPixelX() - pixel->GetPixelX());
+    Int_t distY = TMath::Abs(tmpDig->GetPixelY() - pixel->GetPixelY());
+    if (distX<=1 &&  distY<=1) return 0;
+    if (!iDigit) distance = TMath::Sqrt(distX*distX + distY*distY);
+    else         distance = TMath::Min(distance, TMath::Sqrt(distX*distX + distY*distY));
+  }
+
+  return distance;
+
+}
+
+//====================================================================================================================================================
+
+Bool_t AliMFTCluster::AddPixel(AliMFTDigit *pixel) {
+
+  if (!fIsClusterEditable || fSize>=fNMaxDigitsPerCluster) return kFALSE;
+  if (fSize && (pixel->GetPlane()!=fPlane || pixel->GetDetElemID()!=fDetElemID)) return kFALSE;
+
+  new ((*fDigitsInCluster)[fDigitsInCluster->GetEntries()]) AliMFTDigit(*pixel);
+
+  if (!fSize) {
+    SetPlane(pixel->GetPlane());
+    SetDetElemID(pixel->GetDetElemID());
+  }
+
+  fSize++;
+
+  return kTRUE;
+
+}
+
+//====================================================================================================================================================
+
+void AliMFTCluster::TerminateCluster() {
+
+  Double_t xCenters[fNMaxDigitsPerCluster] = {0};
+  Double_t yCenters[fNMaxDigitsPerCluster] = {0};
+  Double_t nElectrons = 0.;
+
+  for (Int_t iDigit=0; iDigit<fDigitsInCluster->GetEntries(); iDigit++) {
+    AliMFTDigit *tmpDig = (AliMFTDigit*) fDigitsInCluster->At(iDigit);
+    xCenters[iDigit] = tmpDig->GetPixelCenterX();
+    yCenters[iDigit] = tmpDig->GetPixelCenterY();
+    nElectrons      += tmpDig->GetNElectrons();
+    for (Int_t iTrack=0; iTrack<tmpDig->GetNMCTracks(); iTrack++) AddMCLabel(tmpDig->GetMCLabel(iTrack));
+  }
+
+  SetX(TMath::Mean(fDigitsInCluster->GetEntries(), xCenters));
+  SetY(TMath::Mean(fDigitsInCluster->GetEntries(), yCenters));
+  SetZ(((AliMFTDigit*) fDigitsInCluster->At(0))->GetPixelCenterZ());
+
+  Double_t minErrX = ((AliMFTDigit*) fDigitsInCluster->At(0))->GetPixelWidthX() / TMath::Sqrt(12.);
+  Double_t minErrY = ((AliMFTDigit*) fDigitsInCluster->At(0))->GetPixelWidthY() / TMath::Sqrt(12.);
+  Double_t minErrZ = ((AliMFTDigit*) fDigitsInCluster->At(0))->GetPixelWidthZ() / TMath::Sqrt(12.);
+  SetErrX( TMath::Max(TMath::RMS(fDigitsInCluster->GetEntries(), xCenters), minErrX) );
+  SetErrY( TMath::Max(TMath::RMS(fDigitsInCluster->GetEntries(), yCenters), minErrY) );
+  SetErrZ( minErrZ );
+    
+  SetNElectrons(nElectrons);
+
+  fIsClusterEditable = kFALSE;
+
+}
+  
+//====================================================================================================================================================
+
+void AliMFTCluster::AddMCLabel(Int_t label) { 
+
+  if (fNMCTracks>=fNMaxMCTracks) return; 
+
+  for (Int_t iTrack=0; iTrack<fNMCTracks; iTrack++) if (label==fMCLabel[iTrack]) return;
+
+  fMCLabel[fNMCTracks++]=label; 
+
 }
 
 //====================================================================================================================================================
index 898418a..122088b 100644 (file)
 
 #include "AliMUONRawCluster.h"
 #include "AliMUONVCluster.h"
+#include "AliMFTDigit.h"
+#include "TClonesArray.h"
 #include "TObject.h"
+#include "AliMFTConstants.h"
 
 //====================================================================================================================================================
 
@@ -23,26 +26,26 @@ class AliMFTCluster : public TObject {
 public:
 
   AliMFTCluster();
-//   AliMFTCluster(const AliMFTCluster& pt);
-//   AliMFTCluster& operator=(const AliMFTCluster &source);
+  AliMFTCluster(const AliMFTCluster&);
+  AliMFTCluster& operator=(const AliMFTCluster&);
 
   virtual ~AliMFTCluster() {};  // destructor
 
   void SetXYZ(Double_t x, Double_t y, Double_t z) { fX=x; fY=y; fZ=z; }
 
-  void SetX(Double_t x) { fX = x; }
-  void SetY(Double_t y) { fY = y; }
-  void SetZ(Double_t z) { fZ = z; }
+  void SetX(Double_t x) { if(fIsClusterEditable) fX = x; }
+  void SetY(Double_t y) { if(fIsClusterEditable) fY = y; }
+  void SetZ(Double_t z) { if(fIsClusterEditable) fZ = z; }
 
   Double_t GetX() const { return fX; }
   Double_t GetY() const { return fY; }
   Double_t GetZ() const { return fZ; }
   
-  void SetErrXYZ(Double_t errX, Double_t errY, Double_t errZ) { fErrX = errX; fErrY = errY; fErrZ = errZ; }
+  void SetErrXYZ(Double_t errX, Double_t errY, Double_t errZ) { if(fIsClusterEditable) { fErrX = errX; fErrY = errY; fErrZ = errZ; } }
        
-  void SetErrX(Double_t errX) { fErrX = errX; }
-  void SetErrY(Double_t errY) { fErrY = errY; }
-  void SetErrZ(Double_t errZ) { fErrZ = errZ; }
+  void SetErrX(Double_t errX) { if(fIsClusterEditable) fErrX = errX; }
+  void SetErrY(Double_t errY) { if(fIsClusterEditable) fErrY = errY; }
+  void SetErrZ(Double_t errZ) { if(fIsClusterEditable) fErrZ = errZ; }
 
   Double_t GetErrX()  const { return fErrX; }
   Double_t GetErrY()  const { return fErrY; }
@@ -51,17 +54,21 @@ public:
   Double_t GetErrY2() const { return fErrY*fErrY; }
   Double_t GetErrZ2() const { return fErrZ*fErrZ; }
   
-  void     SetNElectrons(Double_t nElectrons) { fNElectrons = nElectrons; }
+  void     SetNElectrons(Double_t nElectrons) { if(fIsClusterEditable) fNElectrons = nElectrons; }
   Double_t GetNElectrons() const { return fNElectrons; }
   
-  void  AddMCLabel(Int_t label) { if (fNMCTracks==fNMaxMCTracks) return; else fMCLabel[fNMCTracks++]=label; }
+  void  AddMCLabel(Int_t label);
   Int_t GetNMCTracks() const { return fNMCTracks; }
   Int_t GetMCLabel(Int_t track) const { if (track<fNMCTracks && track>=0) return fMCLabel[track]; else return -1; }
+  void  SetMCLabel(Int_t track, Int_t labelMC) { if (track<fNMCTracks && track>=0) fMCLabel[track]=labelMC; }
 
-  void  SetPlane(Int_t plane) { fPlane = plane; }
+  void  SetPlane(Int_t plane) { if(fIsClusterEditable) fPlane = plane; }
   Int_t GetPlane() const { return fPlane; }
 
-  void  SetSize(Int_t size) { fSize = size; }
+  void SetDetElemID(Int_t detElemID) { fDetElemID = detElemID; }
+  Int_t GetDetElemID() { return fDetElemID; }
+
+  void  SetSize(Int_t size) { if(fIsClusterEditable) fSize = size; }
   Int_t GetSize() const { return fSize; }
 
   void SetLocalChi2(Double_t chi2) { fLocalChi2 = chi2; }
@@ -70,18 +77,27 @@ public:
   Double_t GetLocalChi2() { return fLocalChi2; }
   Double_t GetTrackChi2() { return fTrackChi2; }
 
+  Bool_t AddPixel(AliMFTDigit *pixel);
+
+  Bool_t IsClusterEditable() { return fIsClusterEditable; }
+  void SetClusterEditable(Bool_t isClusterEditable) { fIsClusterEditable = isClusterEditable; }
+  void TerminateCluster();
+
+  Double_t GetDistanceFromPixel(AliMFTDigit *pixel);
+
   AliMUONRawCluster* CreateMUONCluster();
   
 private:
 
-  static const Int_t fNMaxMCTracks = 30;
-
+  static const Int_t fNMaxMCTracks         = AliMFTConstants::fNMaxMCTracksPerCluster;
+  static const Int_t fNMaxDigitsPerCluster = AliMFTConstants::fNMaxDigitsPerCluster;
+  
   Double_t fX, fY, fZ;   // cluster global coordinates
   Double_t fErrX, fErrY, fErrZ;
 
   Double_t fNElectrons;
   Int_t fNMCTracks;
-  Int_t fPlane;
+  Int_t fPlane, fDetElemID;
   Int_t fMCLabel[fNMaxMCTracks];
 
   Int_t fSize;   // the number of digits composing the cluster
@@ -89,6 +105,10 @@ private:
   Double_t fTrackChi2; // Chi2 of the track when the associated cluster was attached
   Double_t fLocalChi2; // Local chi2 of the associated cluster with respect to the track
   
+  TClonesArray *fDigitsInCluster;   //! (Temporary) Array of the digits composing the cluster
+
+  Bool_t fIsClusterEditable;
+
   ClassDef(AliMFTCluster, 1)
 
 };
index 165ac72..ca22b2e 100644 (file)
 #include "AliMFTSegmentation.h"
 #include "TTree.h"
 #include "TMath.h"
+#include "AliMFTConstants.h"
 #include "AliMFTClusterFinder.h"
 
+const Double_t AliMFTClusterFinder::fCutForAvailableDigits = AliMFTConstants::fCutForAvailableDigits;
+const Double_t AliMFTClusterFinder::fCutForAttachingDigits = AliMFTConstants::fCutForAttachingDigits;
+
+
 ClassImp(AliMFTClusterFinder)
 
 //====================================================================================================================================================
@@ -38,7 +43,8 @@ ClassImp(AliMFTClusterFinder)
 AliMFTClusterFinder::AliMFTClusterFinder() : 
   TObject(),
   fDigitsInCluster(0),
-  fCurrentDig(0),
+  fCurrentDigit(0),
+  fCurrentCluster(0),
   fSegmentation(0),
   fNPlanes(0)
 {
@@ -64,7 +70,7 @@ AliMFTClusterFinder::~AliMFTClusterFinder() {
 
 //====================================================================================================================================================
 
-void AliMFTClusterFinder::Init(const Char_t *nameGeomFile) {
+void AliMFTClusterFinder::Init(Char_t *nameGeomFile) {
 
   fSegmentation = new AliMFTSegmentation(nameGeomFile);
   fNPlanes = fSegmentation -> GetNPlanes();
@@ -97,119 +103,66 @@ void AliMFTClusterFinder::DigitsToClusters(const TObjArray *pDigitList) {
   AliDebug(1, Form("nPlanes = %d", fNPlanes));
 
   StartEvent(); 
+  Bool_t isDigAvailableForNewCluster = kTRUE;
   
   for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
 
     AliDebug(1, Form("Plane %02d", iPlane));
 
     TClonesArray *myDigitList = (TClonesArray*) pDigitList->At(iPlane);
-    //    myDigitList->Sort();
 
     AliDebug(1, Form("myDigitList->GetEntries() = %d", myDigitList->GetEntries()));
 
-    while (myDigitList->GetEntries()) {
-
-      fDigitsInCluster->Clear();
-      
-      Bool_t clusterUpdated=kTRUE;
-
-      //---------------------------------------------------------------------------------------------------------
-
-      while (clusterUpdated) {    // repeat the loop on the digits until no new compatible digit is found
-
-       clusterUpdated = kFALSE;
-
-       for (Int_t iDig=0; iDig<myDigitList->GetEntries(); iDig++) {
-         fCurrentDig = (AliMFTDigit*) myDigitList->At(iDig);
-         if (fDigitsInCluster->GetEntries()<fNMaxDigitsPerCluster) {
-           if (fDigitsInCluster->GetEntries()==0) {
-             new ((*fDigitsInCluster)[fDigitsInCluster->GetEntries()]) AliMFTDigit(*fCurrentDig);
-             myDigitList->Remove(fCurrentDig);
-             clusterUpdated = kTRUE;
-           }
-           else if (IsCurrentDigitCompatible()) {
-             new ((*fDigitsInCluster)[fDigitsInCluster->GetEntries()]) AliMFTDigit(*fCurrentDig);
-             myDigitList->Remove(fCurrentDig);
-             clusterUpdated = kTRUE;
-           }
-         }
-       }
-       
-       if (clusterUpdated) myDigitList->Compress();
-
-      }
-
-      //---------------------------------------------------------------------------------------------------------
-
-      AliDebug(1, "Building new cluster");
-
-      BuildNewCluster(iPlane);  // here the new cluster is built
-
-    }
-
-  }
-
-}
-
-//====================================================================================================================================================
-
-Bool_t AliMFTClusterFinder::IsCurrentDigitCompatible() {
-
-  // where it is decided if the current digit (fCurrentDig) is compatible with the current digits array (fDigitsInCluster)
-
-  for (Int_t iDigit=0; iDigit<fDigitsInCluster->GetEntries(); iDigit++) {
-    AliMFTDigit *tmpDig = (AliMFTDigit*) fDigitsInCluster->At(iDigit);
-    Int_t distX = TMath::Abs(tmpDig->GetPixelX() - fCurrentDig->GetPixelX());
-    Int_t distY = TMath::Abs(tmpDig->GetPixelY() - fCurrentDig->GetPixelY());
-    if (distX<=1 &&  distY<=1) return kTRUE;
-  }
-
-  return kFALSE;
+    Int_t cycleOverDigits = 0;
+    Double_t myCutForAvailableDigits = fCutForAvailableDigits;
     
-}
-
-//====================================================================================================================================================
+    while (myDigitList->GetEntries()) {
 
-void AliMFTClusterFinder::BuildNewCluster(Int_t plane) {
+      for (Int_t iDig=0; iDig<myDigitList->GetEntries(); iDig++) {
 
-  // where a new cluster is built, starting from the array of compatible digits (fDigitsInCluster)
+       AliDebug(1, Form("Check %d: Digit %5d of %5d", cycleOverDigits, iDig, myDigitList->GetEntries()));
 
-  AliDebug(1, Form("Starting cluster building from %d digits", fDigitsInCluster->GetEntries()));
+       fCurrentDigit = (AliMFTDigit*) myDigitList->At(iDig);
+       isDigAvailableForNewCluster = kTRUE;
 
-  AliMFTCluster *newCluster = new AliMFTCluster();
+       for (Int_t iCluster=0; iCluster<fClustersPerPlane[iPlane]->GetEntries(); iCluster++) {
+         fCurrentCluster = (AliMFTCluster*) fClustersPerPlane[iPlane]->At(iCluster);
+         AliDebug(2, Form("Distance between cluster and digit = %f",fCurrentCluster->GetDistanceFromPixel(fCurrentDigit))); 
+         if (fCurrentCluster->GetDistanceFromPixel(fCurrentDigit) < fCutForAttachingDigits) { 
+           fCurrentCluster->AddPixel(fCurrentDigit); 
+           myDigitList->Remove(fCurrentDigit);
+           myDigitList->Compress();
+           iDig--;
+           isDigAvailableForNewCluster = kFALSE;
+           break; 
+         }
+         if (fCurrentCluster->GetDistanceFromPixel(fCurrentDigit) < myCutForAvailableDigits) isDigAvailableForNewCluster=kFALSE;
+       }
 
-  Double_t xCenters[fNMaxDigitsPerCluster] = {0};
-  Double_t yCenters[fNMaxDigitsPerCluster] = {0};
-  Double_t nElectrons = 0.;
+       if (isDigAvailableForNewCluster) {
+         AliMFTCluster *newCluster = new AliMFTCluster();
+         newCluster->AddPixel(fCurrentDigit);
+         myDigitList->Remove(fCurrentDigit);
+         myDigitList->Compress();
+         iDig--;
+         new ((*fClustersPerPlane[iPlane])[fClustersPerPlane[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
+       }
 
-  for (Int_t iDigit=0; iDigit<fDigitsInCluster->GetEntries(); iDigit++) {
-    AliMFTDigit *tmpDig = (AliMFTDigit*) fDigitsInCluster->At(iDigit);
-    xCenters[iDigit] = tmpDig->GetPixelCenterX();
-    yCenters[iDigit] = tmpDig->GetPixelCenterY();
-    nElectrons      += tmpDig->GetNElectrons();
-    for (Int_t iTrack=0; iTrack<tmpDig->GetNMCTracks(); iTrack++) newCluster->AddMCLabel(tmpDig->GetMCLabel(iTrack));
-  }
+      }   // end of cycle over the digits
 
-  newCluster -> SetX(TMath::Mean(fDigitsInCluster->GetEntries(), xCenters));
-  newCluster -> SetY(TMath::Mean(fDigitsInCluster->GetEntries(), yCenters));
-  newCluster -> SetZ(((AliMFTDigit*) fDigitsInCluster->At(0))->GetPixelCenterZ());
+      if (cycleOverDigits) myCutForAvailableDigits -= 0.5;
+      cycleOverDigits++;
 
-  Double_t minErrX = ((AliMFTDigit*) fDigitsInCluster->At(0))->GetPixelWidthX() / TMath::Sqrt(12.);
-  Double_t minErrY = ((AliMFTDigit*) fDigitsInCluster->At(0))->GetPixelWidthY() / TMath::Sqrt(12.);
-  Double_t minErrZ = ((AliMFTDigit*) fDigitsInCluster->At(0))->GetPixelWidthZ() / TMath::Sqrt(12.);
-  newCluster -> SetErrX( TMath::Max(TMath::RMS(fDigitsInCluster->GetEntries(), xCenters), minErrX) );
-  newCluster -> SetErrY( TMath::Max(TMath::RMS(fDigitsInCluster->GetEntries(), yCenters), minErrY) );
-  newCluster -> SetErrZ( minErrZ );
-    
-  newCluster -> SetNElectrons(nElectrons);
-  newCluster -> SetPlane(plane);
+    }   // no more digits to check in current plane!
 
-  newCluster -> SetSize(fDigitsInCluster->GetEntries());
+    for (Int_t iCluster=0; iCluster<fClustersPerPlane[iPlane]->GetEntries(); iCluster++) {
+      fCurrentCluster = (AliMFTCluster*) fClustersPerPlane[iPlane]->At(iCluster);
+      fCurrentCluster -> TerminateCluster();
+    }
 
-  AliDebug(1, Form("Adding cluster #%02d to tree (%f, %f, %f)", 
-                  fClustersPerPlane[plane]->GetEntries(), newCluster->GetX(), newCluster->GetY(), newCluster->GetZ()));
+    AliDebug(1, Form("Found %d clusters in plane %02d", fClustersPerPlane[iPlane]->GetEntries(), iPlane));
 
-  new ((*fClustersPerPlane[plane])[fClustersPerPlane[plane]->GetEntries()]) AliMFTCluster(*newCluster);
+  }  // end of cycle over the planes
 
 }
 
index c464a9f..a8e9b4e 100644 (file)
@@ -20,6 +20,7 @@
 #include "AliMFTSegmentation.h"
 #include "TTree.h"
 #include "TMath.h"
+#include "AliMFTConstants.h"
 
 //====================================================================================================================================================
 
@@ -30,7 +31,7 @@ public:
   AliMFTClusterFinder();
   ~AliMFTClusterFinder();
 
-  void Init(const Char_t *nameGeomFile);
+  void Init(Char_t *nameGeomFile);
   
   void MakeClusterBranch(TTree *treeCluster);
   void SetClusterTreeAddress(TTree *treeCluster);
@@ -39,18 +40,19 @@ public:
   void DigitsToClusters(const TObjArray *pDigitList);
 
   void StartEvent();
-  void BuildNewCluster(Int_t plane);
-  Bool_t IsCurrentDigitCompatible();
 
 private:
  
-  static const Int_t fNMaxDigitsPerCluster = 10;
-  static const Int_t fNMaxPlanes = 20;
+  static const Int_t fNMaxDigitsPerCluster = AliMFTConstants::fNMaxDigitsPerCluster;
+  static const Int_t fNMaxPlanes = AliMFTConstants::fNMaxPlanes;
+  static const Double_t fCutForAvailableDigits;
+  static const Double_t fCutForAttachingDigits;
 
   TClonesArray *fClustersPerPlane[fNMaxPlanes];    // ![fNPlanes] list of clusters [per plane]
 
   TClonesArray *fDigitsInCluster;
-  AliMFTDigit *fCurrentDig;
+  AliMFTDigit *fCurrentDigit;
+  AliMFTCluster *fCurrentCluster;
 
   AliMFTSegmentation *fSegmentation;
 
index f150c38..d93f7b5 100644 (file)
@@ -8,6 +8,7 @@
 #include "AliMFTSegmentation.h"
 #include "TFile.h"
 #include "TH1D.h"
+#include "TH2D.h"
 #include "AliLog.h"
 #include "TString.h"
 
@@ -40,11 +41,12 @@ AliMFTClusterQA::AliMFTClusterQA():
   
   // default constructor
 
-  for (Int_t iPlane=0; iPlane<fMaxNPlanesMFT; iPlane++) {
+  for (Int_t iPlane=0; iPlane<fNMaxPlanes; iPlane++) {
     fHistNClustersPerEvent[iPlane] = 0;
     fHistNPixelsPerCluster[iPlane] = 0;
     fHistClusterSizeX[iPlane] = 0; 
     fHistClusterSizeY[iPlane] = 0;
+    fClusterScatterPlotXY[iPlane] = 0;
   }
 
 }
@@ -87,12 +89,14 @@ Bool_t AliMFTClusterQA::LoadNextEvent() {
   for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) {
     Int_t nClusters = fMFT->GetRecPointsList(iPlane)->GetEntries();
     fHistNClustersPerEvent[iPlane] -> Fill(nClusters);
+    fClusterScatterPlotXY[iPlane]  -> Fill(0., 0.);    // "scaler" bin
     AliDebug(1,Form("nClusters = %5d", nClusters));
     for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
       AliMFTCluster *cluster = (AliMFTCluster*) (fMFT->GetRecPointsList(iPlane))->At(iCluster);
       fHistNPixelsPerCluster[iPlane] -> Fill(cluster->GetSize());
       fHistClusterSizeX[iPlane]      -> Fill(cluster->GetErrX()*1.e4);   // converted in microns
       fHistClusterSizeY[iPlane]      -> Fill(cluster->GetErrY()*1.e4);   // converted in microns
+      fClusterScatterPlotXY[iPlane]  -> Fill(cluster->GetX(), cluster->GetY());
     }
   }
 
@@ -133,6 +137,15 @@ void AliMFTClusterQA::BookHistos() {
     fHistNPixelsPerCluster[iPlane] -> Sumw2();
     fHistClusterSizeX[iPlane]      -> Sumw2();
     fHistClusterSizeY[iPlane]      -> Sumw2();
+
+    //------------------------------------------------------------
+
+    Int_t rMax = Int_t(10.*(fMFT->GetSegmentation()->GetPlane(iPlane)->GetRMaxSupport()));
+    fClusterScatterPlotXY[iPlane] = new TH2D(Form("fClusterScatterPlotXY_Pl%02d",iPlane),
+                                            Form("Cluster scatter plot (Plane%02d)",iPlane),
+                                            2*rMax+1, (-rMax-0.5)/10., (rMax+0.5)/10., 2*rMax+1, (-rMax-0.5)/10., (rMax+0.5)/10.);
+    
+    fClusterScatterPlotXY[iPlane] -> Sumw2();
     
   }
   
@@ -151,6 +164,7 @@ void AliMFTClusterQA::Terminate() {
     fHistNPixelsPerCluster[iPlane] -> Write();
     fHistClusterSizeX[iPlane]      -> Write();
     fHistClusterSizeY[iPlane]      -> Write();
+    fClusterScatterPlotXY[iPlane]  -> Write();
   }
 
   fFileOut -> Close();
index 0301e23..b1cb3a5 100644 (file)
@@ -11,6 +11,7 @@
 #include "AliMFTSegmentation.h"
 #include "TFile.h"
 #include "TH1D.h"
+#include "TH2D.h"
 #include "AliLog.h"
 #include "TString.h"
 
@@ -41,12 +42,11 @@ private:
 
 protected:
 
-  static const Int_t fMaxNPlanesMFT = 20;
+  static const Int_t fNMaxPlanes = AliMFTConstants::fNMaxPlanes;
 
-  TH1D *fHistNClustersPerEvent[fMaxNPlanesMFT], *fHistNPixelsPerCluster[fMaxNPlanesMFT];
-  TH1D *fHistClusterSizeX[fMaxNPlanesMFT], *fHistClusterSizeY[fMaxNPlanesMFT];
-
-  TClonesArray *fMFTClusterArray[fMaxNPlanesMFT];
+  TH1D *fHistNClustersPerEvent[fNMaxPlanes], *fHistNPixelsPerCluster[fNMaxPlanes];
+  TH1D *fHistClusterSizeX[fNMaxPlanes], *fHistClusterSizeY[fNMaxPlanes];
+  TH2D *fClusterScatterPlotXY[fNMaxPlanes];
 
   AliLoader *fMFTLoader;
   AliRunLoader *fRunLoader;
diff --git a/MFT/AliMFTConstants.cxx b/MFT/AliMFTConstants.cxx
new file mode 100644 (file)
index 0000000..d2c3217
--- /dev/null
@@ -0,0 +1,45 @@
+ /**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+//====================================================================================================================================================
+//
+//      Constants for the Muon Forward Tracker
+//
+//      Contact author: antonio.uras@cern.ch
+//
+//====================================================================================================================================================
+
+#include "TClass.h"
+#include "AliMFTConstants.h"
+
+ClassImp(AliMFTConstants)
+
+const Double_t AliMFTConstants::fCutForAvailableDigits = 5.;
+const Double_t AliMFTConstants::fCutForAttachingDigits = 1.;
+
+const Double_t AliMFTConstants::fElossPerElectron = 3.62e-09;
+
+const Double_t AliMFTConstants::fRadiusMin = 2.225;
+
+const Double_t AliMFTConstants::fActiveSuperposition = 0.05;
+                                 
+const Double_t AliMFTConstants::fHeightActive = 0.5;
+const Double_t AliMFTConstants::fHeightReadout = 0.3;
+
+const Double_t AliMFTConstants::fSupportExtMargin = fHeightReadout + 0.3;
+
+const Double_t AliMFTConstants::fRadLengthSi = 9.37;
+
+//====================================================================================================================================================
diff --git a/MFT/AliMFTConstants.h b/MFT/AliMFTConstants.h
new file mode 100644 (file)
index 0000000..44ceb27
--- /dev/null
@@ -0,0 +1,58 @@
+#ifndef AliMFTConstants_H
+#define AliMFTConstants_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//====================================================================================================================================================
+//
+//      Constants for the Muon Forward Tracker
+//
+//      Contact author: antonio.uras@cern.ch
+//
+//====================================================================================================================================================
+
+#include <TObject.h>
+
+class AliMFTConstants : public TObject {
+
+public:
+  
+  static const Int_t fNMaxPlanes = 20; 
+
+  static const Int_t fNMaxDigitsPerCluster = 12;       ///< max number of digits per cluster
+  static const Double_t fCutForAvailableDigits;   ///<
+  static const Double_t fCutForAttachingDigits;   ///<
+
+  static const Int_t fNMaxMCTracksPerCluster = 30;   ///< max number of MC tracks sharing the same MFT cluster
+  static const Int_t fNMaxMCTracksPerDigit = 10;     ///< max number of MC tracks sharing the same MFT digit
+
+  static const Double_t fElossPerElectron;
+
+  // minimum radial distance of the MFT sensors. To be carefully coordinated with fActiveSuperposition
+  static const Double_t fRadiusMin;  ///<
+        
+  // superposition between the active elements tasselling the MFT planes, for having a full acceptance coverage even in case of 10 degrees inclined tracks
+  static const Double_t fActiveSuperposition;  ///<
+                                                
+  static const Double_t fHeightActive;   ///< height of the active elements
+  static const Double_t fHeightReadout;  ///< height of the readout elements attached to the active ones
+        
+  // minimum border size between the end of the support plane and the sensors: fHeightReadout + 0.3
+  static const Double_t fSupportExtMargin;  ///<
+
+  static const Int_t fNMaxDetElemPerPlane = 1000;  ///<
+
+  static const Double_t fRadLengthSi;    ///< expressed in cm
+
+protected:
+
+  AliMFTConstants() : TObject() {}
+  virtual ~AliMFTConstants(){}
+
+  ClassDef(AliMFTConstants, 0)   // MFT global constants 
+
+};
+       
+#endif
+
index 9364fc1..3390cb1 100644 (file)
 //====================================================================================================================================================
 
 #include "AliDigit.h"
+#include "AliMFTConstants.h"
 #include "AliMFTDigit.h"
 
+const Double_t AliMFTDigit::fElossPerElectron = AliMFTConstants::fElossPerElectron;
+
 ClassImp(AliMFTDigit)
 
 //====================================================================================================================================================
@@ -48,8 +51,18 @@ AliMFTDigit::AliMFTDigit():
 
   // default cosntructor
 
-  for (Int_t iTr=0; iTr<fNMaxMCTracks; iTr++) fMCLabel[iTr] = -1; 
+  for (Int_t iTrack=0; iTrack<fNMaxMCTracksPerDigit; iTrack++) fMCLabel[iTrack] = -1;
+
+}
+
+//====================================================================================================================================================
+
+void  AliMFTDigit::AddMCLabel(Int_t label) { 
+
+  if (fNMCTracks<0 || fNMCTracks>=fNMaxMCTracksPerDigit) return; 
+  fMCLabel[fNMCTracks++] = label;
 
 }
 
 //====================================================================================================================================================
+
index 3dd5c19..715e45a 100644 (file)
@@ -13,6 +13,7 @@
 //====================================================================================================================================================
 
 #include "AliDigit.h"
+#include "AliMFTConstants.h"
 
 //====================================================================================================================================================
 
@@ -21,7 +22,6 @@ class AliMFTDigit: public AliDigit {
 public:
 
   AliMFTDigit();
-//   AliMFTDigit(const AliMFTDigit &d);
 
   virtual ~AliMFTDigit() {}
 
@@ -42,9 +42,10 @@ public:
   }
   void SetEloss(Double_t sig) { fEloss = sig; fNElectrons = fEloss/fElossPerElectron; }
 
-  void  AddMCLabel(Int_t label) { if (fNMCTracks>=fNMaxMCTracks) return; else fMCLabel[fNMCTracks++]=label; }
+  void  AddMCLabel(Int_t label); 
+
   Int_t GetNMCTracks() const { return fNMCTracks; }
-  Int_t GetMCLabel(Int_t track) const { if (track<fNMCTracks && track>=0) return fMCLabel[track]; else return -1; }
+  Int_t GetMCLabel(Int_t track) const { if (track<fNMCTracks && track>=0 && fNMCTracks>0) return fMCLabel[track]; else return -1; }
 
   Double_t  GetEloss()         const { return fEloss; }
   Double_t  GetNElectrons()    const { return fNElectrons; }
@@ -62,9 +63,9 @@ public:
     
 protected:
     
-  static const Double_t fElossPerElectron = 3.62e-09;
-  static const Int_t fNMaxMCTracks = 10;
-
+  static const Double_t fElossPerElectron;
+  static const Int_t fNMaxMCTracksPerDigit = AliMFTConstants::fNMaxMCTracksPerDigit;
+  
   Int_t fNMCTracks;
 
   Int_t fPixelX;
@@ -81,7 +82,7 @@ protected:
   Double_t fEloss;       // total signal as Eloss in the medium
   Double_t fNElectrons; 
   
-  Int_t fMCLabel[fNMaxMCTracks];
+  Int_t fMCLabel[fNMaxMCTracksPerDigit];
 
   ClassDef(AliMFTDigit,3)
 
@@ -90,6 +91,3 @@ protected:
 //====================================================================================================================================================
 
 #endif
-
-
-
index f6ba7c1..8bd6dbf 100644 (file)
@@ -50,7 +50,7 @@ AliMFTDigitizer::AliMFTDigitizer():
 
 //====================================================================================================================================================
 
-AliMFTDigitizer::AliMFTDigitizer(AliDigitizationInput *digInp) :
+AliMFTDigitizer::AliMFTDigitizer(AliDigitizationInput *digInp):
   AliDigitizer(digInp),
   fNPlanes(0),
   fSegmentation(0)
@@ -65,7 +65,7 @@ void AliMFTDigitizer::Digitize(Option_t*) {
   // This method is responsible for merging sdigits to a list of digits
 
   AliDebug(1, "************************************************************************");
-  AliDebug(1, "************************ AliMFTDigitizer::Exec *************************");
+  AliDebug(1, "************************ AliMFTDigitizer::Digitize *********************");
   AliDebug(1, "************************************************************************");
   
   if (!fSegmentation) {
index e1272b9..1189299 100644 (file)
@@ -41,8 +41,8 @@ public:
   
 protected:
  
-  static const Int_t fNMaxPlanes = 20;     // max number of MFT planes
-  static const Int_t fNMaxMCTracks = 10;   // max MC tracks sharing a digit  
+  static const Int_t fNMaxPlanes   = AliMFTConstants::fNMaxPlanes;             // max number of MFT planes
+  static const Int_t fNMaxMCTracks = AliMFTConstants::fNMaxMCTracksPerDigit;   // max MC tracks sharing a digit  
 
   Int_t fNPlanes;    
   
index 45f81b6..7b001dc 100644 (file)
Binary files a/MFT/AliMFTGeometry.root and b/MFT/AliMFTGeometry.root differ
index 78fb56f..0e9cc07 100644 (file)
@@ -24,7 +24,6 @@
 #include "TNamed.h"
 #include "THnSparse.h"
 #include "TClonesArray.h"
-#include "AliMFTPlane.h"
 #include "TAxis.h"
 #include "TPave.h"
 #include "TCanvas.h"
 #include "TEllipse.h"
 #include "TMath.h"
 #include "AliLog.h"
+#include "AliMFTConstants.h"
+#include "AliMFTPlane.h"
+
+const Double_t AliMFTPlane::fRadiusMin           = AliMFTConstants::fRadiusMin;
+const Double_t AliMFTPlane::fActiveSuperposition = AliMFTConstants::fActiveSuperposition;
+const Double_t AliMFTPlane::fHeightActive        = AliMFTConstants::fHeightActive;
+const Double_t AliMFTPlane::fHeightReadout       = AliMFTConstants::fHeightReadout;
+const Double_t AliMFTPlane::fSupportExtMargin    = AliMFTConstants::fSupportExtMargin;
 
 ClassImp(AliMFTPlane)
 
@@ -39,7 +46,7 @@ ClassImp(AliMFTPlane)
 
 AliMFTPlane::AliMFTPlane():
   TNamed(),
-  fPlaneNumber(0),
+  fPlaneNumber(-1),
   fZCenter(0), 
   fRMinSupport(0), 
   fRMax(0),
@@ -59,12 +66,6 @@ AliMFTPlane::AliMFTPlane():
   fSupportElements(0)
 {
 
-  fPlaneNumber = -1;
-
-  fActiveElements  = new TClonesArray("THnSparseC");
-  fReadoutElements = new TClonesArray("THnSparseC");
-  fSupportElements = new TClonesArray("THnSparseC");
-
   // default constructor
 
 }
@@ -73,7 +74,7 @@ AliMFTPlane::AliMFTPlane():
 
 AliMFTPlane::AliMFTPlane(const Char_t *name, const Char_t *title):
   TNamed(name, title),
-  fPlaneNumber(0),
+  fPlaneNumber(-1),
   fZCenter(0), 
   fRMinSupport(0), 
   fRMax(0),
@@ -88,18 +89,12 @@ AliMFTPlane::AliMFTPlane(const Char_t *name, const Char_t *title):
   fEquivalentSilicon(0),
   fEquivalentSiliconBeforeFront(0),
   fEquivalentSiliconBeforeBack(0),
-  fActiveElements(0),
-  fReadoutElements(0),
-  fSupportElements(0)
+  fActiveElements(new TClonesArray("THnSparseC")),
+  fReadoutElements(new TClonesArray("THnSparseC")),
+  fSupportElements(new TClonesArray("THnSparseC"))
 {
 
-  fPlaneNumber = -1;
-
-  fActiveElements  = new TClonesArray("THnSparseC");
-  fReadoutElements = new TClonesArray("THnSparseC");
-  fSupportElements = new TClonesArray("THnSparseC");
-
-  // default constructor
+  // constructor
 
 }
 
@@ -122,12 +117,16 @@ AliMFTPlane::AliMFTPlane(const AliMFTPlane& plane):
   fEquivalentSilicon(plane.fEquivalentSilicon),
   fEquivalentSiliconBeforeFront(plane.fEquivalentSiliconBeforeFront),
   fEquivalentSiliconBeforeBack(plane.fEquivalentSiliconBeforeBack),
-  fActiveElements(plane.fActiveElements),
-  fReadoutElements(plane.fReadoutElements),
-  fSupportElements(plane.fSupportElements)
+  fActiveElements(new TClonesArray("THnSparseC")),
+  fReadoutElements(new TClonesArray("THnSparseC")),
+  fSupportElements(new TClonesArray("THnSparseC"))
 {
 
   // copy constructor
+
+  *fActiveElements  = *plane.fActiveElements;
+  *fReadoutElements = *plane.fReadoutElements;
+  *fSupportElements = *plane.fSupportElements;
   
 }
 
@@ -138,35 +137,33 @@ AliMFTPlane& AliMFTPlane::operator=(const AliMFTPlane& plane) {
   // Assignment operator
 
   // check assignement to self
-  if (this == &plane) return *this;
+  if (this != &plane) {
 
-  // base class assignement
-  TNamed::operator=(plane);
-  
-  // clear memory
-  Clear();
-
-  fPlaneNumber                  = plane.fPlaneNumber;
-  fZCenter                      = plane.fZCenter; 
-  fRMinSupport                  = plane.fRMinSupport; 
-  fRMax                         = plane.fRMax;
-  fRMaxSupport                  = plane.fRMaxSupport;
-  fPixelSizeX                   = plane.fPixelSizeX;
-  fPixelSizeY                   = plane.fPixelSizeY; 
-  fThicknessActive              = plane.fThicknessActive; 
-  fThicknessSupport             = plane.fThicknessSupport; 
-  fThicknessReadout             = plane.fThicknessReadout;
-  fZCenterActiveFront           = plane.fZCenterActiveFront;
-  fZCenterActiveBack            = plane.fZCenterActiveBack;
-  fEquivalentSilicon            = plane.fEquivalentSilicon;
-  fEquivalentSiliconBeforeFront = plane.fEquivalentSiliconBeforeFront;
-  fEquivalentSiliconBeforeBack  = plane.fEquivalentSiliconBeforeBack;
-  fActiveElements               = plane.fActiveElements;
-  fReadoutElements              = plane.fReadoutElements;
-  fSupportElements              = plane.fSupportElements;
+    // base class assignement
+    TNamed::operator=(plane);
+    
+    fPlaneNumber                  = plane.fPlaneNumber;
+    fZCenter                      = plane.fZCenter; 
+    fRMinSupport                  = plane.fRMinSupport; 
+    fRMax                         = plane.fRMax;
+    fRMaxSupport                  = plane.fRMaxSupport;
+    fPixelSizeX                   = plane.fPixelSizeX;
+    fPixelSizeY                   = plane.fPixelSizeY; 
+    fThicknessActive              = plane.fThicknessActive; 
+    fThicknessSupport             = plane.fThicknessSupport; 
+    fThicknessReadout             = plane.fThicknessReadout;
+    fZCenterActiveFront           = plane.fZCenterActiveFront;
+    fZCenterActiveBack            = plane.fZCenterActiveBack;
+    fEquivalentSilicon            = plane.fEquivalentSilicon;
+    fEquivalentSiliconBeforeFront = plane.fEquivalentSiliconBeforeFront;
+    fEquivalentSiliconBeforeBack  = plane.fEquivalentSiliconBeforeBack;
+    *fActiveElements              = *plane.fActiveElements;
+    *fReadoutElements             = *plane.fReadoutElements;
+    *fSupportElements             = *plane.fSupportElements;
+  }
 
   return *this;
-
+  
 }
 
 //====================================================================================================================================================
@@ -206,7 +203,7 @@ Bool_t AliMFTPlane::Init(Int_t    planeNumber,
   fRMax = fRMinSupport + (fHeightActive-fActiveSuperposition) * 
     (Int_t((fRMax-fRMinSupport-fHeightActive)/(fHeightActive-fActiveSuperposition))+1) + fHeightActive;
   
-  fRMaxSupport = TMath::Sqrt(fHeightActive*(rMax-fHeightActive) + fRMax*fRMax) + fSupportExtMargin;
+  fRMaxSupport = TMath::Sqrt(fHeightActive*(2.*rMax-fHeightActive) + fRMax*fRMax) + fSupportExtMargin;
  
   return kTRUE;
  
@@ -397,7 +394,7 @@ THnSparseC* AliMFTPlane::GetSupportElement(Int_t id) {
 
 //====================================================================================================================================================
 
-void AliMFTPlane::DrawPlane(Char_t *opt) {
+void AliMFTPlane::DrawPlane(Option_t *opt) {
 
   // ------------------- "FRONT" option ------------------
 
@@ -415,7 +412,7 @@ void AliMFTPlane::DrawPlane(Char_t *opt) {
     h->SetYTitle("y [cm]");
     h->Draw();
 
-    printf("Created hist\n");
+    AliInfo("Created hist");
 
     TEllipse *supportExt = new TEllipse(0.0, 0.0, fRMaxSupport, fRMaxSupport);
     TEllipse *supportInt = new TEllipse(0.0, 0.0, fRMinSupport, fRMinSupport);
index 6961ee2..4ae87b0 100644 (file)
@@ -48,9 +48,9 @@ public:
   
   Bool_t CreateStructure();
 
-  Int_t GetNActiveElements()  { return fActiveElements->GetEntries();  }
-  Int_t GetNReadoutElements() { return fReadoutElements->GetEntries(); }
-  Int_t GetNSupportElements() { return fSupportElements->GetEntries(); }
+  Int_t GetNActiveElements()  const { return fActiveElements->GetEntries();  }
+  Int_t GetNReadoutElements() const { return fReadoutElements->GetEntries(); }
+  Int_t GetNSupportElements() const { return fSupportElements->GetEntries(); }
 
   TClonesArray* GetActiveElements()  { return fActiveElements;  }
   TClonesArray* GetReadoutElements() { return fReadoutElements; }
@@ -60,37 +60,37 @@ public:
   THnSparseC* GetReadoutElement(Int_t id);
   THnSparseC* GetSupportElement(Int_t id);
 
-  Bool_t IsFront(THnSparseC *element) { return (element->GetAxis(2)->GetXmin() < fZCenter); }
+  Bool_t IsFront(THnSparseC *element) const { return (element->GetAxis(2)->GetXmin() < fZCenter); }
 
-  void DrawPlane(Char_t *opt="");
+  void DrawPlane(Option_t *opt="");
 
-  Double_t GetRMinSupport() { return fRMinSupport; }
-  Double_t GetRMaxSupport() { return fRMaxSupport; }
+  Double_t GetRMinSupport() const { return fRMinSupport; }
+  Double_t GetRMaxSupport() const { return fRMaxSupport; }
   Double_t GetThicknessSupport() { return GetSupportElement(0)->GetAxis(2)->GetXmax() - GetSupportElement(0)->GetAxis(2)->GetXmin(); }
   
-  Double_t GetZCenter() { return fZCenter; }
-  Double_t GetZCenterActiveFront() { return fZCenterActiveFront; }
-  Double_t GetZCenterActiveBack()  { return fZCenterActiveBack; }
+  Double_t GetZCenter()            const { return fZCenter; }
+  Double_t GetZCenterActiveFront() const { return fZCenterActiveFront; }
+  Double_t GetZCenterActiveBack()  const { return fZCenterActiveBack; }
 
   void SetEquivalentSilicon(Double_t equivalentSilicon)                       { fEquivalentSilicon            = equivalentSilicon; }
   void SetEquivalentSiliconBeforeFront(Double_t equivalentSiliconBeforeFront) { fEquivalentSiliconBeforeFront = equivalentSiliconBeforeFront; }
   void SetEquivalentSiliconBeforeBack(Double_t equivalentSiliconBeforeBack)   { fEquivalentSiliconBeforeBack  = equivalentSiliconBeforeBack; }
-  Double_t GetEquivalentSilicon()            { return fEquivalentSilicon; }
-  Double_t GetEquivalentSiliconBeforeFront() { return fEquivalentSiliconBeforeFront; }
-  Double_t GetEquivalentSiliconBeforeBack()  { return fEquivalentSiliconBeforeBack; }
+  Double_t GetEquivalentSilicon()            const { return fEquivalentSilicon; }
+  Double_t GetEquivalentSiliconBeforeFront() const { return fEquivalentSiliconBeforeFront; }
+  Double_t GetEquivalentSiliconBeforeBack()  const { return fEquivalentSiliconBeforeBack; }
   
 private:
 
   // measures in cm
 
-  static const Double_t fRadiusMin = 2.225;           // minimum radial distance of the MFT sensors. To be carefully coordinated with fActiveSuperposition
+  static const Double_t fRadiusMin;            // minimum radial distance of the MFT sensors. To be carefully coordinated with fActiveSuperposition
 
-  static const Double_t fActiveSuperposition = 0.05;  // superposition between the active elements tasselling the MFT planes, for having a 
-                                                      // full acceptance coverage even in case of 10 degrees inclined tracks
-  static const Double_t fHeightActive = 0.5;          // height of the active elements
-  static const Double_t fHeightReadout = 0.3;         // height of the readout elements attached to the active ones
+  static const Double_t fActiveSuperposition;  // superposition between the active elements tasselling the MFT planes, for having a 
+                                               // full acceptance coverage even in case of 10 degrees inclined tracks
+  static const Double_t fHeightActive;         // height of the active elements
+  static const Double_t fHeightReadout;        // height of the readout elements attached to the active ones
 
-  static const Double_t fSupportExtMargin = 0.3 + 0.3;    // minimum border size between the end of the support plane and the sensors: fHeightReadout + 0.3
+  static const Double_t fSupportExtMargin;     // minimum border size between the end of the support plane and the sensors: fHeightReadout + 0.3
 
   Int_t fPlaneNumber;
 
index b8198b6..8d955df 100644 (file)
@@ -66,7 +66,7 @@ void AliMFTReconstructor::Init() {
 
   fDigits = new TObjArray(fNPlanes);
   fDigits->SetOwner(kTRUE);
-  for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) fDigits->AddAt(new TClonesArray("AliMFTDigit",fNMaxDigitPerPlane),iPlane);
+  for (Int_t iPlane=0; iPlane<fNPlanes; iPlane++) fDigits->AddAt(new TClonesArray("AliMFTDigit"),iPlane);
 
   AliInfo("    ************* Using the MFT reconstructor! ****** ");
 
index 9299765..d0ec953 100644 (file)
@@ -43,8 +43,6 @@ private:
   AliMFTReconstructor(const AliMFTReconstructor&);              // Not implemented
   AliMFTReconstructor &operator=(const AliMFTReconstructor&);   // Not implemented
 
-  static const Int_t fNMaxDigitPerPlane = 10000;
-
   TObjArray  *fDigits;     
   Int_t      fNPlanes;
 
index deb7e22..fc89866 100644 (file)
@@ -37,22 +37,18 @@ AliMFTSegmentation::AliMFTSegmentation():
   fMFTPlanes(0)
 { 
 
-  // TO BE CHECKED
-  
   // default constructor
 
-  fMFTPlanes = new TClonesArray("AliMFTPlane", fNMaxPlanes);
-
 }
 
 //====================================================================================================================================================
 
 AliMFTSegmentation::AliMFTSegmentation(const Char_t *nameGeomFile): 
   TObject(),
-  fMFTPlanes(0)
+  fMFTPlanes(new TClonesArray("AliMFTPlane", fNMaxPlanes))
 { 
 
-  fMFTPlanes = new TClonesArray("AliMFTPlane", fNMaxPlanes);
+  // constructor
 
   Float_t zCenter, rMin, rMax, pixelSizeX, pixelSizeY, thicknessActive, thicknessSupport, thicknessReadout;
   Float_t equivalentSilicon, equivalentSiliconBeforeFront, equivalentSiliconBeforeBack;
@@ -102,7 +98,7 @@ AliMFTSegmentation::AliMFTSegmentation(const Char_t *nameGeomFile):
 
 //====================================================================================================================================================
 
-THnSparseC* AliMFTSegmentation::GetDetElem(Int_t detElemID) {
+THnSparseC* AliMFTSegmentation::GetDetElem(Int_t detElemID) const {
       
   // Find det elem
 
index c5c7d5f..9600a3d 100644 (file)
@@ -18,6 +18,7 @@
 #include "TMath.h"
 #include "AliMFTPlane.h"
 #include "TMath.h"
+#include "AliMFTConstants.h"
 
 //====================================================================================================================================================
 
@@ -30,34 +31,34 @@ public:
 
   virtual ~AliMFTSegmentation() {}
 
-  THnSparseC* GetDetElem(Int_t detElemID);
+  THnSparseC* GetDetElem(Int_t detElemID) const;
 
-  Int_t GetDetElemID(Int_t plane, Int_t detElem) { return fNMaxDetElemPerPlane*plane + detElem; }
+  Int_t GetDetElemID(Int_t plane, Int_t detElem) const { return fNMaxDetElemPerPlane*plane + detElem; }
     
   Bool_t Hit2PixelID(Double_t xHit, Double_t yHit, Int_t detElemID, Int_t &xPixel, Int_t &yPixel);  
 
-  Double_t GetPixelSizeX(Int_t detElemID) { THnSparseC *detElem = GetDetElem(detElemID); return detElem->GetAxis(0)->GetBinWidth(1); }
-  Double_t GetPixelSizeY(Int_t detElemID) { THnSparseC *detElem = GetDetElem(detElemID); return detElem->GetAxis(1)->GetBinWidth(1); }
-  Double_t GetPixelSizeZ(Int_t detElemID) { THnSparseC *detElem = GetDetElem(detElemID); return detElem->GetAxis(2)->GetBinWidth(1); }
+  Double_t GetPixelSizeX(Int_t detElemID) const { THnSparseC *detElem = GetDetElem(detElemID); return detElem->GetAxis(0)->GetBinWidth(1); }
+  Double_t GetPixelSizeY(Int_t detElemID) const { THnSparseC *detElem = GetDetElem(detElemID); return detElem->GetAxis(1)->GetBinWidth(1); }
+  Double_t GetPixelSizeZ(Int_t detElemID) const { THnSparseC *detElem = GetDetElem(detElemID); return detElem->GetAxis(2)->GetBinWidth(1); }
 
-  Double_t GetPixelCenterX(Int_t detElemID, Int_t iPixel) { THnSparseC *detElem = GetDetElem(detElemID); return detElem->GetAxis(0)->GetBinCenter(iPixel+1); }
-  Double_t GetPixelCenterY(Int_t detElemID, Int_t iPixel) { THnSparseC *detElem = GetDetElem(detElemID); return detElem->GetAxis(1)->GetBinCenter(iPixel+1); }
-  Double_t GetPixelCenterZ(Int_t detElemID, Int_t iPixel) { THnSparseC *detElem = GetDetElem(detElemID); return -1.*(detElem->GetAxis(2)->GetBinCenter(iPixel+1)); }
+  Double_t GetPixelCenterX(Int_t detElemID, Int_t iPixel) const { THnSparseC *detElem = GetDetElem(detElemID); return detElem->GetAxis(0)->GetBinCenter(iPixel+1); }
+  Double_t GetPixelCenterY(Int_t detElemID, Int_t iPixel) const { THnSparseC *detElem = GetDetElem(detElemID); return detElem->GetAxis(1)->GetBinCenter(iPixel+1); }
+  Double_t GetPixelCenterZ(Int_t detElemID, Int_t iPixel) const { THnSparseC *detElem = GetDetElem(detElemID); return -1.*(detElem->GetAxis(2)->GetBinCenter(iPixel+1)); }
 
-  Int_t GetNPlanes() { return fMFTPlanes->GetEntries(); }
+  Int_t GetNPlanes() const { return fMFTPlanes->GetEntries(); }
 
-  AliMFTPlane* GetPlane(Int_t iPlane) { if (iPlane>=0 && iPlane<fMFTPlanes->GetEntries()) return (AliMFTPlane*) fMFTPlanes->At(iPlane); else return NULL; }
+  AliMFTPlane* GetPlane(Int_t iPlane) const { if (iPlane>=0 && iPlane<fMFTPlanes->GetEntries()) return (AliMFTPlane*) fMFTPlanes->At(iPlane); else return NULL; }
  
 protected:
 
-  static const Int_t fNMaxPlanes = 20;                // max number of MFT planes
-  static const Double_t fRadiusMin = 2.225;           // minimum radial distance of the MFT sensors. To be carefully coordinated with fDetElemSuperposition
-  static const Double_t fDetElemSuperposition = 0.05; // superposition between bands tasselling the MFT planes, for having a full acceptance coverage
-                                                      // even in case of 10 degrees inclined tracks
-  static const Double_t fHeightDetElem = 0.5;         // height of the active volume bands composing the planes
-  static const Double_t fSupportExtMargin = 0.3;      // minimum border size between the end of the support plane and the sensors
+  static const Int_t fNMaxPlanes = AliMFTConstants::fNMaxPlanes;                // max number of MFT planes
+  static const Double_t fRadiusMin;             // minimum radial distance of the MFT sensors. To be carefully coordinated with fDetElemSuperposition
+  static const Double_t fDetElemSuperposition;  // superposition between bands tasselling the MFT planes, for having a full acceptance coverage
+                                                // even in case of 10 degrees inclined tracks
+  static const Double_t fHeightDetElem;         // height of the active volume bands composing the planes
+  static const Double_t fSupportExtMargin;      // minimum border size between the end of the support plane and the sensors
 
-  static const Int_t fNMaxDetElemPerPlane = 1000;
+  static const Int_t fNMaxDetElemPerPlane = AliMFTConstants::fNMaxDetElemPerPlane;
 
   TClonesArray *fMFTPlanes;
 
index 0a5897f..3e2a6a6 100644 (file)
@@ -31,6 +31,7 @@
 #include "TMatrixD.h"
 #include "TParticle.h"
 #include "AliMuonForwardTrack.h"
+#include "AliMFTConstants.h"
 
 ClassImp(AliMuonForwardTrack)
 
@@ -40,11 +41,17 @@ AliMuonForwardTrack::AliMuonForwardTrack():
   AliMUONTrack(),
   fMUONTrack(0),
   fMCTrackRef(0),
-  fMFTClusters(0)
+  fMFTClusters(0),
+  fNWrongClustersMC(-1)
 {
 
   // default constructor
-  for (Int_t iPlane=0; iPlane<fMaxNPlanesMFT; iPlane++) fPlaneExists[iPlane] = kFALSE;
+  
+  for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = kFALSE;
+  for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
+    fParentMCLabel[iParent] = -1;
+    fParentPDGCode[iParent] =  0;
+  }
   fMFTClusters = new TClonesArray("AliMFTCluster");
 
 }
@@ -55,11 +62,16 @@ AliMuonForwardTrack::AliMuonForwardTrack(AliMUONTrack *MUONTrack):
   AliMUONTrack(),
   fMUONTrack(0),
   fMCTrackRef(0),
-  fMFTClusters(0)
+  fMFTClusters(0),
+  fNWrongClustersMC(-1)
 {
 
   SetMUONTrack(MUONTrack);
-  for (Int_t iPlane=0; iPlane<fMaxNPlanesMFT; iPlane++) fPlaneExists[iPlane] = kFALSE;
+  for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = kFALSE;
+  for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
+    fParentMCLabel[iParent] = -1;
+    fParentPDGCode[iParent] =  0;
+  }
   fMFTClusters = new TClonesArray("AliMFTCluster");
 
 }
@@ -68,13 +80,21 @@ AliMuonForwardTrack::AliMuonForwardTrack(AliMUONTrack *MUONTrack):
 
 AliMuonForwardTrack::AliMuonForwardTrack(const AliMuonForwardTrack& track): 
   AliMUONTrack(track),
-  fMUONTrack(track.fMUONTrack),
-  fMCTrackRef(track.fMCTrackRef),
-  fMFTClusters(track.fMFTClusters)
+  fMUONTrack(0x0),
+  fMCTrackRef(0x0),
+  fMFTClusters(0x0),
+  fNWrongClustersMC(track.fNWrongClustersMC)
 {
 
   // copy constructor
-  for (Int_t iPlane=0; iPlane<fMaxNPlanesMFT; iPlane++) fPlaneExists[iPlane] = (track.fPlaneExists)[iPlane];
+  fMUONTrack        = new AliMUONTrack(*(track.fMUONTrack));
+  if (track.fMCTrackRef) fMCTrackRef = new TParticle(*(track.fMCTrackRef));
+  fMFTClusters      = new TClonesArray(*(track.fMFTClusters));
+  for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = (track.fPlaneExists)[iPlane];
+  for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
+    fParentMCLabel[iParent] = (track.fParentMCLabel)[iParent];
+    fParentPDGCode[iParent] = (track.fParentPDGCode)[iParent];
+  }
   
 }
 
@@ -93,11 +113,16 @@ AliMuonForwardTrack& AliMuonForwardTrack::operator=(const AliMuonForwardTrack& t
   // clear memory
   Clear();
   
-  fMUONTrack    = track.fMUONTrack;
-  fMCTrackRef   = track.fMCTrackRef;
-  fMFTClusters  = track.fMFTClusters;
-
-  for (Int_t iPlane=0; iPlane<fMaxNPlanesMFT; iPlane++) fPlaneExists[iPlane] = (track.fPlaneExists)[iPlane];
+  fMUONTrack        = new AliMUONTrack(*(track.fMUONTrack));
+  if (track.fMCTrackRef) fMCTrackRef = new TParticle(*(track.fMCTrackRef));
+  fMFTClusters      = new TClonesArray(*(track.fMFTClusters));
+  fNWrongClustersMC = track.fNWrongClustersMC;
+
+  for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) fPlaneExists[iPlane] = (track.fPlaneExists)[iPlane];
+  for (Int_t iParent=0; iParent<fgkNParentsMax; iParent++) {
+    fParentMCLabel[iParent] = (track.fParentMCLabel)[iParent];
+    fParentPDGCode[iParent] = (track.fParentPDGCode)[iParent];
+  }
   
   return *this;
 
@@ -113,6 +138,7 @@ void AliMuonForwardTrack::SetMUONTrack(AliMUONTrack *MUONTrack) {
   }
 
   fMUONTrack = MUONTrack;
+  SetGlobalChi2(fMUONTrack->GetGlobalChi2());
   
 }
 
@@ -133,8 +159,8 @@ void AliMuonForwardTrack::SetMCTrackRef(TParticle *MCTrackRef) {
 
 void AliMuonForwardTrack::AddTrackParamAtMFTCluster(AliMUONTrackParam &trackParam, AliMFTCluster &mftCluster) {
 
-  AliDebug(1, Form("fMFTClusters=%p has %d entries", fMFTClusters, fMFTClusters->GetEntries()));
-  Int_t iMFTCluster = fMFTClusters->GetEntries();
+  AliDebug(1, Form("Before adding: this->fMFTClusters=%p has %d entries", this->fMFTClusters, this->fMFTClusters->GetEntries()));
+  Int_t iMFTCluster = this->fMFTClusters->GetEntries();
   AliDebug(1, Form("mftCluster->GetX() = %f  mftCluster->GetY() = %f  mftCluster->GetErrX() = %f  cmftCluster->GetErrY() = %f", 
           mftCluster.GetX(), mftCluster.GetY(), mftCluster.GetErrX(), mftCluster.GetErrY()));
   AliMUONVCluster *muonCluster = (AliMUONVCluster*) mftCluster.CreateMUONCluster();
@@ -142,15 +168,18 @@ void AliMuonForwardTrack::AddTrackParamAtMFTCluster(AliMUONTrackParam &trackPara
   trackParam.SetUniqueID(iMFTCluster);    // we profit of this slot to store the reference to the corresponding MFTCluster
   AliDebug(1, Form("Now adding muonCluster %p and trackParam %p",muonCluster, &trackParam));
   AddTrackParamAtCluster(trackParam, *muonCluster, kTRUE);
-  AliDebug(1, Form("GetTrackParamAtCluster() = %p  has %d entries",GetTrackParamAtCluster(), GetTrackParamAtCluster()->GetEntries()));
   // we pass the parameters this->GetTrackParamAtCluster()->First() to the Kalman Filter algorithm: they will be updated!!
   Double_t chi2Kalman = RunKalmanFilter(*(AliMUONTrackParam*)(GetTrackParamAtCluster()->First()));
+  AliDebug(1, Form("Adding Kalman chi2 = %f to global chi2 = %f", chi2Kalman, GetGlobalChi2()));
   Double_t newGlobalChi2 = GetGlobalChi2() + chi2Kalman;
   mftCluster.SetLocalChi2(chi2Kalman);
   mftCluster.SetTrackChi2(newGlobalChi2);
-  new ((*fMFTClusters)[iMFTCluster]) AliMFTCluster(mftCluster);
+  new ((*(this->fMFTClusters))[iMFTCluster]) AliMFTCluster(mftCluster);
+  AliDebug(1, Form("GetTrackParamAtCluster() = %p  has %d entries while this->fMFTClusters=%p has %d entries",
+                  GetTrackParamAtCluster(), GetTrackParamAtCluster()->GetEntries(), this->fMFTClusters, this->fMFTClusters->GetEntries()));
   AliDebug(1, Form("muonCluster->GetZ() = %f, trackParam->GetZ() = %f",muonCluster->GetZ(), trackParam.GetZ()));
   SetGlobalChi2(newGlobalChi2);
+  AliDebug(1, Form("New global chi2 = %f", GetGlobalChi2()));
   ((AliMUONTrackParam*) GetTrackParamAtCluster()->First())->SetTrackChi2(newGlobalChi2);
   for (Int_t iPar=0; iPar<GetTrackParamAtCluster()->GetEntries(); iPar++) {
     AliDebug(1, Form("GetTrackParamAtCluster()->At(%d)->GetClusterPtr() = %p",
@@ -210,12 +239,12 @@ AliMUONVCluster* AliMuonForwardTrack::GetMUONCluster(Int_t iMUONCluster) {
 AliMFTCluster* AliMuonForwardTrack::GetMFTCluster(Int_t iMFTCluster) {
 
   if (iMFTCluster<0 || iMFTCluster>=GetNMFTClusters()) {
-    AliError("Invalid MFT cluster index. NULL pointer will be returned");
+    AliError(Form("Invalid MFT cluster index (%d). GetNMFTClusters()=%d. NULL pointer will be returned", iMFTCluster, GetNMFTClusters()));
     return NULL;
   }
 
   AliMUONTrackParam *trackParam = GetTrackParamAtMFTCluster(iMFTCluster);
-  AliMFTCluster *mftCluster = (AliMFTCluster*) fMFTClusters->At(trackParam->GetUniqueID());
+  AliMFTCluster *mftCluster = (AliMFTCluster*) this->fMFTClusters->At(trackParam->GetUniqueID());
 
   return mftCluster;
 
index 7546db3..0b9af18 100644 (file)
@@ -20,6 +20,7 @@
 #include "TMatrixD.h"
 #include "TClonesArray.h"
 #include "TParticle.h"
+#include "AliMFTConstants.h"
 
 //====================================================================================================================================================
 
@@ -27,6 +28,8 @@ class AliMuonForwardTrack : public AliMUONTrack {
 
 public:
 
+  static const Int_t fgkNParentsMax =  5;   ///< maximum number of parents
+
   AliMuonForwardTrack();
   AliMuonForwardTrack(AliMUONTrack *MUONTrack);
 
@@ -50,7 +53,7 @@ public:
   Bool_t PlaneExists(Int_t iPlane) { return fPlaneExists[iPlane]; }
 
   Int_t GetNMUONClusters() { return fMUONTrack->GetNClusters(); }
-  Int_t GetNMFTClusters()  { return GetNClusters(); }
+  Int_t GetNMFTClusters()  { return fMFTClusters->GetEntries(); }
 
   Int_t GetMCLabelMUONTrack() { return fMUONTrack->GetMCLabel(); }
 
@@ -63,17 +66,33 @@ public:
   Double_t GetOffsetX(Double_t x, Double_t z);
   Double_t GetOffsetY(Double_t y, Double_t z);
 
+  void SetParentMCLabel(Int_t iParent, Int_t MClabel) { if (0<=iParent && iParent<fgkNParentsMax) fParentMCLabel[iParent] = MClabel; }
+  void SetParentPDGCode(Int_t iParent, Int_t PDGCode) { if (0<=iParent && iParent<fgkNParentsMax) fParentPDGCode[iParent] = PDGCode; }
+
+  Int_t GetParentMCLabel(Int_t iParent) { if (0<=iParent && iParent<fgkNParentsMax) return fParentMCLabel[iParent]; else return -1; }
+  Int_t GetParentPDGCode(Int_t iParent) { if (0<=iParent && iParent<fgkNParentsMax) return fParentPDGCode[iParent]; else return  0; }
+
+  void SetNWrongClustersMC(Int_t nClusters) { fNWrongClustersMC = nClusters; }
+  Int_t GetNWrongClustersMC() { return fNWrongClustersMC; }
+
+  Double_t Pt() { return TMath::Sqrt(TMath::Power(GetTrackParamAtMFTCluster(0)->Px(),2)+TMath::Power(GetTrackParamAtMFTCluster(0)->Py(),2)); }
+  
 protected:
 
-  static const Int_t fMaxNPlanesMFT = 20;
+  static const Int_t fNMaxPlanes = AliMFTConstants::fNMaxPlanes;        // max number of MFT planes
 
-  Bool_t fPlaneExists[fMaxNPlanesMFT];
+  Bool_t fPlaneExists[fNMaxPlanes];
 
   AliMUONTrack *fMUONTrack;
   TParticle *fMCTrackRef;
 
   TClonesArray *fMFTClusters;
 
+  Int_t fParentMCLabel[fgkNParentsMax];    ///< MC label of parents and grandparents
+  Int_t fParentPDGCode[fgkNParentsMax];    ///< PDG code of parents and grandparents 
+
+  Int_t fNWrongClustersMC;    // number of wrong associated MC clusters
+
   ClassDef(AliMuonForwardTrack,1)
     
 };
index 8a4db65..0ee7433 100644 (file)
@@ -1,12 +1,15 @@
 //================================================================================================================================
 
 void AliMuonForwardTrackFinder(Int_t run=0,
+                              Int_t matching=0,
                               const Char_t *readDir= ".",
                               const Char_t *outDir = ".",
                               Int_t nEventsToAnalyze = -1) {
   
   //  TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1, AliMagF::k5kG));
 
+  // AliLog::SetClassDebugLevel("AliMuonForwardTrackFinder", 1);
+
   AliCDBManager* man = AliCDBManager::Instance();
   man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
   man->SetSpecificStorage("GRP/GRP/Data", Form("local://%s",gSystem->pwd()));
@@ -24,7 +27,7 @@ void AliMuonForwardTrackFinder(Int_t run=0,
   //  finder -> SetLowPtCut(0.5);
   finder -> SetExtrapOriginTransvError(0.05);   // To be imposed if the gen. vertex is fixed in (0,0,0)
   finder -> SetGaussianBlurZVert(5.0);          // To be imposed if the gen. vertex is fixed in (0,0,0) 
-  finder -> SetMatchingMode(0);                 // 0 -> real matching   1 -> ideal matching
+  finder -> SetMatchingMode(matching);          // 0 -> real matching   1 -> ideal matching
   finder -> SetMinResearchRadiusAtLastPlane(0.5);
 
   while (finder->LoadNextTrack()) continue;
index 66fb665..ffa44e7 100644 (file)
@@ -47,6 +47,7 @@
 #include "AliMFTCluster.h"
 #include "AliMFT.h"
 #include "AliMFTSegmentation.h"
+#include "AliMFTConstants.h"
 
 #include "AliMuonForwardTrackFinder.h"
 
@@ -59,6 +60,8 @@
 //
 //====================================================================================================================================================
 
+const Double_t AliMuonForwardTrackFinder::fRadLengthSi = AliMFTConstants::fRadLengthSi;
+
 ClassImp(AliMuonForwardTrackFinder)
 
 //=====================================================================================================
@@ -97,16 +100,15 @@ AliMuonForwardTrackFinder::AliMuonForwardTrackFinder():
   fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane(0),
   fHistDistanceGoodClusterFromTrackAtLastPlane(0),
 
-  fNtuFinalCandidates1(0),
-  fNtuFinalBestCandidates1(0),
-  fNtuFinalCandidates2(0),
-  fNtuFinalBestCandidates2(0),
+  fNtuFinalCandidates(0),
+  fNtuFinalBestCandidates(0),
 
   fCanvas(0),
 
   fTxtMuonHistory(0), 
   fTxtTrackGoodClusters(0), 
-  fTxtTrackFinalChi2(0), 
+  fTxtTrackFinalChi2(0),
+  fTxtTrackMomentum(0),
   fTxtFinalCandidates(0), 
   fTxtDummy(0),
   fTxtAllClust(0), 
@@ -118,7 +120,8 @@ AliMuonForwardTrackFinder::AliMuonForwardTrackFinder():
   fMrkClustMC(0), 
   fMrkClustOfTrack(0),
 
-  fCountRealTracksAnalyzed(0), 
+  fCountRealTracksAnalyzed(0),
+  fMaxNTracksToBeAnalyzed(99999999),
   fCountRealTracksWithRefMC(0), 
   fCountRealTracksWithRefMC_andTrigger(0),
   fCountRealTracksWithRefMC_andTrigger_andGoodPt(0),
@@ -136,6 +139,7 @@ AliMuonForwardTrackFinder::AliMuonForwardTrackFinder():
   fMFTClusterTree(0),
   fMuonTrackReco(0),
   fCurrentTrack(0),
+  fFinalBestCandidate(0),
   fIsCurrentMuonTrackable(0),
   fCandidateTracks(0),
   fTrackStore(0),
@@ -145,6 +149,7 @@ AliMuonForwardTrackFinder::AliMuonForwardTrackFinder():
   fMFT(0),
   fSegmentation(0),
   fOutputTreeFile(0),
+  fOutputQAFile(0),
   fOutputEventTree(0),
   fMuonForwardTracks(0),
   fMatchingMode(-1),
@@ -156,7 +161,7 @@ AliMuonForwardTrackFinder::AliMuonForwardTrackFinder():
 
   // Default constructor
 
-  for (Int_t iPlane=0; iPlane<fMaxNPlanesMFT; iPlane++) {
+  for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) {
 
     fHistNTracksAfterExtrapolation[iPlane] = 0;
     fHistChi2Cluster_GoodCluster[iPlane] = 0;
@@ -168,14 +173,14 @@ AliMuonForwardTrackFinder::AliMuonForwardTrackFinder():
     fHistChi2Cluster_GoodCluster[iPlane] = 0;
     fHistChi2Cluster_BadCluster[iPlane]  = 0;
     
-    fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] = 0;
-    fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] = 0;
+    fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] = 0;
+    fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] = 0;
     
     fZPlane[iPlane] = 0.;
     fRPlaneMax[iPlane] = 0.;
     fRPlaneMin[iPlane] = 0.;
     
-    for (Int_t i=0; i<4; i++) fGrMFTPlane[iPlane][i] = 0;
+    for (Int_t i=0; i<4; i++) fGrMFTPlane[i][iPlane] = 0;
     fCircleExt[iPlane] = 0;
     fCircleInt[iPlane] = 0;
     
@@ -199,7 +204,7 @@ AliMuonForwardTrackFinder::AliMuonForwardTrackFinder():
   fMFTClusterTree = 0;
   fCandidateTracks = 0;
 
-  fOutputTreeFile = new TFile("MuonGlobalTracks.root", "recreate");
+  fOutputTreeFile    = new TFile("MuonGlobalTracks.root", "recreate");
   fOutputEventTree   = new TTree("AliMuonForwardTracks", "Tree of AliMuonForwardTracks");
   fMuonForwardTracks = new TClonesArray("AliMuonForwardTrack");
   fOutputEventTree   -> Branch("tracks", &fMuonForwardTracks);
@@ -208,25 +213,120 @@ AliMuonForwardTrackFinder::AliMuonForwardTrackFinder():
 
 //=====================================================================================================
 
+AliMuonForwardTrackFinder::~AliMuonForwardTrackFinder() {
+
+  for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) {
+
+    delete fHistNTracksAfterExtrapolation[iPlane];
+    delete fHistChi2Cluster_GoodCluster[iPlane];
+    delete fHistChi2Cluster_BadCluster[iPlane];
+    delete fHistResearchRadius[iPlane];
+    
+    delete fHistChi2Cluster_GoodCluster[iPlane];
+    delete fHistChi2Cluster_BadCluster[iPlane];
+    
+    delete fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane];
+    delete fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane];
+    
+    for (Int_t i=0; i<4; i++) delete fGrMFTPlane[i][iPlane];
+    delete fCircleExt[iPlane];
+    delete fCircleInt[iPlane];
+    
+    delete fTxtTrackChi2[iPlane];
+    
+    delete fMFTClusterArray[iPlane];
+    delete fMFTClusterArrayFront[iPlane];
+    delete fMFTClusterArrayBack[iPlane];
+
+  }
+
+  delete fNtuFinalCandidates;
+  delete fNtuFinalBestCandidates;
+
+  delete fHistPtSpectrometer;
+  delete fHistPtMuonTrackWithGoodMatch;
+  delete fHistPtMuonTrackWithBadMatch;
+  delete fHistRadiusEndOfAbsorber;
+
+  delete fHistNGoodClustersForFinalTracks; 
+  delete fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane;                //
+  delete fHistDistanceGoodClusterFromTrackAtLastPlane;                                         //
+
+  delete fCanvas;
+
+  delete fTxtMuonHistory;
+  delete fTxtTrackGoodClusters;
+  delete fTxtTrackFinalChi2;
+  delete fTxtTrackMomentum;
+  delete fTxtFinalCandidates;
+  delete fTxtDummy;
+  delete fTxtAllClust;
+  delete fTxtClustGoodChi2;
+  delete fTxtClustMC;
+  delete fTxtClustOfTrack;
+  delete fMrkAllClust;
+  delete fMrkClustGoodChi2;
+  delete fMrkClustMC;
+  delete fMrkClustOfTrack;
+  delete fFileCluster;
+  delete fFileESD;
+  delete fFile_gAlice;
+
+  delete fRunLoader;
+  delete fMFTLoader;
+  delete fMuonRecoCheck;
+
+  delete fMFTClusterTree;
+
+  delete fMuonTrackReco;
+  delete fCurrentTrack;
+  delete fFinalBestCandidate;
+
+  delete fCandidateTracks;
+
+  delete fTrackStore;
+  delete fTrackRefStore;
+  
+  delete fNextTrack;
+  
+  delete fStack;
+
+  delete fMFT;
+  delete fSegmentation;
+
+  delete fOutputTreeFile; 
+  delete fOutputQAFile;
+  delete fOutputEventTree;
+
+  delete fMuonForwardTracks;
+
+  delete fGRPData;
+  delete fRunInfo;       
+
+}
+
+//=====================================================================================================
+
 void AliMuonForwardTrackFinder::Init(Int_t nRun, 
                                     Char_t *readDir,
                                     Char_t *outDir,
                                     Int_t nEventsToAnalyze) {
   
   if (fRunLoader) {
-    printf("WARNING: run already initialized!!\n");
+    AliInfo("WARNING: run already initialized!!\n");
   }
 
   SetRun(nRun);
   SetReadDir(readDir);
   SetOutDir(outDir);
 
-  printf("input  dir = %s\n", fReadDir.Data());
-  printf("output dir = %s\n", fOutDir.Data());
+  AliInfo(Form("input  dir = %s\n", fReadDir.Data()));
+  AliInfo(Form("output dir = %s\n", fOutDir.Data()));
 
   // -------------------------- initializing files...
 
-  printf("initializing files for run %d...\n", fRun);
+  AliInfo(Form("initializing files for run %d...\n", fRun));
 
   Char_t geoFileName[300];
   Char_t esdFileName[300];
@@ -242,19 +342,19 @@ void AliMuonForwardTrackFinder::Init(Int_t nRun,
   if (!gGeoManager) {
     TGeoManager::Import(geoFileName);
     if (!gGeoManager) {
-      printf("getting geometry from file %s failed", geoFileName);
+      AliError(Form("getting geometry from file %s failed", geoFileName));
       return;
     }
   }
   
   fFileESD = new TFile(esdFileName);
   if (!fFileESD || !fFileESD->IsOpen()) return;
-  else printf("file %s successfully opened\n", fFileESD->GetName());
+  else AliInfo(Form("file %s successfully opened\n", fFileESD->GetName()));
   
   fMuonRecoCheck = new AliMUONRecoCheck(esdFileName, Form("%s/generated/", fReadDir.Data()));       // Utility class to check reconstruction
   fFile_gAlice = new TFile(gAliceName);
   if (!fFile_gAlice || !fFile_gAlice->IsOpen()) return;
-  else printf("file %s successfully opened\n", fFile_gAlice->GetName());
+  else AliInfo(Form("file %s successfully opened\n", fFile_gAlice->GetName()));
   
   fRunLoader = AliRunLoader::Open(gAliceName);
   gAlice = fRunLoader->GetAliRun();
@@ -279,29 +379,28 @@ void AliMuonForwardTrackFinder::Init(Int_t nRun,
   
   fMFTClusterTree = fMFTLoader->TreeR();
 
-
   Int_t nEventsInFile = fMuonRecoCheck->NumberOfEvents();
   if (!nEventsInFile) {
-    printf("no events available!!!\n");
+    AliError("no events available!!!\n");
     return;
   }
   if (nEventsInFile<nEventsToAnalyze || nEventsToAnalyze<0) fNEventsToAnalyze = nEventsInFile;
   else fNEventsToAnalyze = nEventsToAnalyze;
 
-  fCandidateTracks = new TClonesArray("AliMuonForwardTrack",200000);
+  fCandidateTracks = new TClonesArray("AliMuonForwardTrack",50000);
 
   // -------------------------- initializing histograms...
 
-  printf("\ninitializing histograms...\n");
+  AliInfo("\ninitializing histograms...\n");
   BookHistos();
   SetTitleHistos();
-  printf("... done!\n\n");
+  AliInfo("... done!\n\n");
 
   // -------------------------- initializing graphics...
 
-  printf("initializing graphics...\n");
+  AliInfo("initializing graphics...\n");
   BookPlanes();
-  printf("... done!\n\n");
+  AliInfo("... done!\n\n");
 
   SetSigmaSpectrometerCut(4.0);
   SetSigmaClusterCut(4.5);
@@ -325,7 +424,7 @@ Bool_t AliMuonForwardTrackFinder::LoadNextEvent() {
 
   fCountRealTracksAnalyzedOfEvent = 0;
   
-  printf(" **** analyzing event # %d  \n", fEv);
+  AliInfo(Form(" **** analyzing event # %d  \n", fEv));
   
   fTrackStore    = fMuonRecoCheck->ReconstructedTracks(fEv);
   fTrackRefStore = fMuonRecoCheck->ReconstructibleTracks(fEv);
@@ -333,7 +432,7 @@ Bool_t AliMuonForwardTrackFinder::LoadNextEvent() {
   fRunLoader->GetEvent(fEv);
   if (!fMFTLoader->TreeR()->GetEvent()) return kFALSE;
   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
-    printf("plane %02d: nClusters = %d\n", iPlane, (fMFT->GetRecPointsList(iPlane))->GetEntries());
+    AliDebug(1, Form("plane %02d: nClusters = %d\n", iPlane, (fMFT->GetRecPointsList(iPlane))->GetEntries()));
     fMFTClusterArray[iPlane] = fMFT->GetRecPointsList(iPlane);
   }
   SeparateFrontBackClusters();
@@ -357,13 +456,19 @@ Int_t AliMuonForwardTrackFinder::LoadNextTrack() {
 
   // load next muon track from the reconstructed event
 
+  if (fCountRealTracksAnalyzed>fMaxNTracksToBeAnalyzed) return kFALSE;
+
   if (!fCountRealTracksAnalyzed) if (!LoadNextEvent()) return kFALSE;
+  if (fCountRealTracksAnalyzed==fMaxNTracksToBeAnalyzed) {
+    fCountRealTracksAnalyzed++;
+    if (!LoadNextEvent()) return kFALSE;
+  }
 
   while ( !(fMuonTrackReco = static_cast<AliMUONTrack*>(fNextTrack->Next())) ) if (!LoadNextEvent()) return kFALSE;
 
-  printf("**************************************************************************************\n");
-  printf("***************************   MUON TRACK %3d   ***************************************\n", fCountRealTracksAnalyzedOfEvent);
-  printf("**************************************************************************************\n");
+  AliDebug(1, "**************************************************************************************\n");
+  AliDebug(1, Form("***************************   MUON TRACK %3d   ***************************************\n", fCountRealTracksAnalyzedOfEvent));
+  AliDebug(1, "**************************************************************************************\n");
 
   fCountRealTracksAnalyzed++;
 
@@ -398,14 +503,12 @@ Int_t AliMuonForwardTrackFinder::LoadNextTrack() {
 
   CheckCurrentMuonTrackable();
 
-  PrintParticleHistory();
-  
   if (fMuonTrackReco->GetMatchTrigger()) fCountRealTracksWithRefMC_andTrigger++;
   
   // the track we are going to build, starting from fMuonTrackReco and adding the MFT clusters
   AliMuonForwardTrack *track = new ((*fCandidateTracks)[0]) AliMuonForwardTrack();
   track -> SetMUONTrack(fMuonTrackReco);
-  if (fLabelMC>=0) track -> SetMCTrackRef(fStack->Particle(fLabelMC));
+  if (fLabelMC>=0 && fStack->Particle(fLabelMC)) track -> SetMCTrackRef(fStack->Particle(fLabelMC));
   track -> SetMCLabel(fMuonTrackReco->GetMCLabel());
   track -> SetMatchTrigger(fMuonTrackReco->GetMatchTrigger());
   
@@ -415,6 +518,7 @@ Int_t AliMuonForwardTrackFinder::LoadNextTrack() {
   Double_t thetaSpectrometer = TMath::ATan(ptSpectrometer/param->Pz());
   if (thetaSpectrometer<0.) thetaSpectrometer += TMath::Pi();
   Double_t etaSpectrometer = -1.*TMath::Log(TMath::Tan(0.5*thetaSpectrometer));
+  //  fOutputQAFile->cd();
   fHistPtSpectrometer -> Fill(ptSpectrometer);
   
   // if the transverse momentum in the Muon Spectrometer is smaller than the threshold, skip to the next track
@@ -426,6 +530,7 @@ Int_t AliMuonForwardTrackFinder::LoadNextTrack() {
   Double_t xEndOfAbsorber = trackParamEndOfAbsorber.GetNonBendingCoor();
   Double_t yEndOfAbsorber = trackParamEndOfAbsorber.GetBendingCoor();
   Double_t rAbsorber      = TMath::Sqrt(xEndOfAbsorber*xEndOfAbsorber + yEndOfAbsorber*yEndOfAbsorber);
+  //  fOutputQAFile->cd();
   fHistRadiusEndOfAbsorber -> Fill(rAbsorber);
   
   // if the radial distance of the track at the end of the absorber is smaller than a radius corresponding to 
@@ -454,13 +559,16 @@ Int_t AliMuonForwardTrackFinder::LoadNextTrack() {
        }
       }
       fCandidateTracks->Compress();
-      if (fIsCurrentMuonTrackable) fHistNTracksAfterExtrapolation[iPlane] -> Fill(fCandidateTracks->GetEntriesFast());
+      if (fIsCurrentMuonTrackable) {
+       //      fOutputQAFile->cd();
+       fHistNTracksAfterExtrapolation[iPlane] -> Fill(fCandidateTracks->GetEntriesFast());
+      }
     }
 
-    else if (fMatchingMode==kIdealMatching) {
+    else if (fMatchingMode==kIdealMatching && fIsCurrentMuonTrackable) {
       fCurrentTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(0);
-      printf("plane %02d: fCandidateTracks->GetEntriesFast() = %d   fCandidateTracks->UncheckedAt(0) = %p   fCurrentTrack = %p\n", 
-            iPlane, fCandidateTracks->GetEntriesFast(), fCandidateTracks->UncheckedAt(0), fCurrentTrack);
+      AliDebug(2, Form("plane %02d: fCandidateTracks->GetEntriesFast() = %d   fCandidateTracks->UncheckedAt(0) = %p   fCurrentTrack = %p\n", 
+                      iPlane, fCandidateTracks->GetEntriesFast(), fCandidateTracks->UncheckedAt(0), fCurrentTrack));
       AttachGoodClusterInPlane(iPlane);
     }
 
@@ -468,19 +576,53 @@ Int_t AliMuonForwardTrackFinder::LoadNextTrack() {
   
   // -------------------------- END OF THE CYCLE OVER THE MFT PLANES --------------------------------------------
   
+  AliDebug(1, "Finished cycle over planes");
+
+  Double_t momentum = ptSpectrometer * TMath::CosH(etaSpectrometer);
+  fTxtTrackMomentum = new TLatex(0.10, 0.70, Form("P_{spectro} = %3.1f GeV/c", momentum));
+
   if (fMatchingMode==kIdealMatching) {
-    printf("Adding track to output tree...\n");
+    AliDebug(1, "Adding track to output tree...\n");
+    fFinalBestCandidate = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(0);
     AliMuonForwardTrack *newTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(0);
-    new ((*fMuonForwardTracks)[fMuonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack);  // AU
-    printf("...track added!\n");
+    new ((*fMuonForwardTracks)[fMuonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack);
+    AliDebug(1, "...track added!\n");
     fCandidateTracks->Clear();
     fCountRealTracksAnalyzedOfEvent++;
+    fCountRealTracksAnalyzedWithFinalCandidates++;
+    PrintParticleHistory();
+    FillPlanesWithTrackHistory();
+
+    Double_t chi2AtPlane[fNMaxPlanes] = {0};
+    Int_t nGoodClusters = 0;
+    Int_t nMFTClusters  = fFinalBestCandidate->GetNMFTClusters();
+//     Int_t nMUONClusters = fFinalBestCandidate->GetNMUONClusters();
+    Int_t plane = 0;
+    for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) {
+      while (!fFinalBestCandidate->PlaneExists(plane)) plane++;
+      AliMFTCluster *localCluster = fFinalBestCandidate->GetMFTCluster(iCluster);
+      chi2AtPlane[plane] = localCluster->GetLocalChi2();
+      if (IsCorrectMatch(localCluster)) nGoodClusters++;
+//       Int_t nClustersGlobalTrack = nMUONClusters + (nMFTClusters-iCluster);        // Muon Spectrometer clusters + clusters in the Vertex Telescope
+//       Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
+//       chi2AtPlane[plane] /= Double_t(ndfGlobalTrack);
+      plane++;
+    }
+    for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
+      fTxtTrackChi2[iPlane] = new TLatex(0.55*fRPlaneMax[fNPlanesMFT-1], 
+                                        0.90*fRPlaneMax[fNPlanesMFT-1], 
+                                        Form("#chi^{2} = %3.1f", chi2AtPlane[iPlane]));
+    }
+    fTxtTrackFinalChi2 = new TLatex(0.20, 0.44, Form("#chi^{2}_{final} = %3.1f", chi2AtPlane[0]));
+
+    if (fDrawOption) DrawPlanes();
     return 5;
   }
-
+  
   // If we have several final tracks, we must find the best candidate:
 
   Int_t nFinalTracks = fCandidateTracks->GetEntriesFast();
+  AliDebug(1, Form("nFinalTracks = %d", nFinalTracks));
 
   if (nFinalTracks) fCountRealTracksAnalyzedWithFinalCandidates++;
   
@@ -488,8 +630,12 @@ Int_t AliMuonForwardTrackFinder::LoadNextTrack() {
   Bool_t bestCandidateExists       = kFALSE;
   Int_t nGoodClustersBestCandidate = 0;
   Int_t idBestCandidate            = 0;
-  Double_t chi2HistoryForBestCandidate[fMaxNPlanesMFT]={0};    // chi2 on each plane, for the best candidate
-  for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) chi2HistoryForBestCandidate[iPlane] = -1.;
+  Double_t chi2HistoryForBestCandidate[fNMaxPlanes] = {0};  // chi2 on each plane, for the best candidate
+  Double_t nClustersPerPlane[fNMaxPlanes] = {0};
+  for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
+    chi2HistoryForBestCandidate[iPlane] = -1.;
+    nClustersPerPlane[iPlane] = fMFTClusterArray[iPlane] -> GetEntries();
+  }
   
   fTxtFinalCandidates = new TLatex(0.10, 0.78, Form("N_{FinalCandidates} = %d", nFinalTracks));
   
@@ -500,61 +646,64 @@ Int_t AliMuonForwardTrackFinder::LoadNextTrack() {
     
     AliMuonForwardTrack *finalTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iTrack);
     
-    Double_t chi2AtPlane[fMaxNPlanesMFT]={0};
+    Double_t chi2AtPlane[fNMaxPlanes] = {0};
     Int_t nGoodClusters = 0;
     Int_t nMFTClusters  = finalTrack->GetNMFTClusters();
-    Int_t nMUONClusters = finalTrack->GetNMUONClusters();
+//     Int_t nMUONClusters = finalTrack->GetNMUONClusters();
 
     Int_t plane = 0;
     for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) {
       while (!finalTrack->PlaneExists(plane)) plane++;
       AliMFTCluster *localCluster = finalTrack->GetMFTCluster(iCluster);
-      chi2AtPlane[plane++] = localCluster->GetTrackChi2();
+      chi2AtPlane[plane] = localCluster->GetLocalChi2();
       if (IsCorrectMatch(localCluster)) nGoodClusters++;
-      Int_t nClustersGlobalTrack = nMUONClusters + (nMFTClusters-iCluster);        // Muon Spectrometer clusters + clusters in the Vertex Telescope
-      Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
-      chi2AtPlane[plane] /= Double_t(ndfGlobalTrack);
+//       Int_t nClustersGlobalTrack = nMUONClusters + (nMFTClusters-iCluster);        // Muon Spectrometer clusters + clusters in the Vertex Telescope
+//       Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
+//       chi2AtPlane[plane] /= Double_t(ndfGlobalTrack);
+      plane++;
     }
     
-    if (fIsCurrentMuonTrackable) fHistNGoodClustersForFinalTracks -> Fill(nGoodClusters);
-
-    fNtuFinalCandidates1 -> Fill(Double_t(fRun),
-                                Double_t(fEv),
-                                Double_t(fCountRealTracksAnalyzedOfEvent),
-                                Double_t(nFinalTracks),
-                                Double_t(nClustersMC),
-                                Double_t(nGoodClusters),
-                                ptSpectrometer,
-                                thetaSpectrometer,
-                                etaSpectrometer);
-
-    fNtuFinalCandidates2 -> Fill(chi2AtPlane[0],
-                                chi2AtPlane[1],
-                                chi2AtPlane[2],
-                                chi2AtPlane[3],
-                                chi2AtPlane[4],
-                                chi2AtPlane[5],
-                                chi2AtPlane[6],
-                                chi2AtPlane[7],
-                                chi2AtPlane[8],
-                                chi2AtPlane[9],
-                                chi2AtPlane[10],
-                                chi2AtPlane[11],
-                                chi2AtPlane[12],
-                                chi2AtPlane[13],
-                                chi2AtPlane[14]);
+    if (fIsCurrentMuonTrackable) {
+      //      fOutputQAFile->cd();
+      fHistNGoodClustersForFinalTracks -> Fill(nGoodClusters);
+    }
+
+    //    fOutputQAFile->cd();
+
+    Float_t finalCandidatesInfo[] = {Double_t(fRun),
+                                    Double_t(fEv),
+                                    Double_t(fCountRealTracksAnalyzedOfEvent),
+                                    Double_t(nFinalTracks),
+                                    Double_t(nClustersMC),
+                                    Double_t(nGoodClusters),
+                                    ptSpectrometer,
+                                    thetaSpectrometer,
+                                    etaSpectrometer, 
+                                    chi2AtPlane[0],
+                                    chi2AtPlane[1],
+                                    chi2AtPlane[2],
+                                    chi2AtPlane[3],
+                                    chi2AtPlane[4],
+                                    chi2AtPlane[5],
+                                    chi2AtPlane[6],
+                                    chi2AtPlane[7],
+                                    chi2AtPlane[8]};
     
+    fNtuFinalCandidates -> Fill(finalCandidatesInfo);
+
     // now comparing the tracks with various criteria, in order to find the best one
     
     Double_t theVariable = 0.;
+//     theVariable = chi2AtPlane[0];
     for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) theVariable += chi2AtPlane[iCluster];
     theVariable /= Double_t(nMFTClusters);
+    
       
     if (theVariable<theVariable_Best || theVariable_Best<0.) {
       nGoodClustersBestCandidate = nGoodClusters;
       for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) chi2HistoryForBestCandidate[iPlane] = chi2AtPlane[iPlane];
       theVariable_Best = theVariable;
-      fTxtTrackFinalChi2 = new TLatex(0.20, 0.52, Form("#chi^{2}_{final} = %3.1f", chi2HistoryForBestCandidate[0]));
+      fTxtTrackFinalChi2 = new TLatex(0.20, 0.44, Form("#chi^{2}_{final} = %3.1f", chi2HistoryForBestCandidate[0]));
       for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
        fTxtTrackChi2[iPlane] = new TLatex(0.55*fRPlaneMax[fNPlanesMFT-1], 
                                           0.90*fRPlaneMax[fNPlanesMFT-1], 
@@ -567,45 +716,56 @@ Int_t AliMuonForwardTrackFinder::LoadNextTrack() {
     // ----------------------------------------------------------
 
   }
-  
+
   if (nFinalTracks) {
-    FillPlanesWithTrackHistory((AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(idBestCandidate));
+    fFinalBestCandidate = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(idBestCandidate);
+    AliInfo(Form("fFinalBestCandidate->GetNMFTClusters() = %d\n",  fFinalBestCandidate->GetNMFTClusters()));
+    PrintParticleHistory();
+    FillPlanesWithTrackHistory();
     AliMuonForwardTrack *newTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(idBestCandidate);
-    new ((*fMuonForwardTracks)[fMuonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack);   // AU
-    fNtuFinalBestCandidates1 -> Fill(Double_t(fRun),
-                                    Double_t(fEv),
-                                    Double_t(fCountRealTracksAnalyzedOfEvent),
-                                    Double_t(nFinalTracks),
-                                    Double_t(nClustersMC),
-                                    Double_t(nGoodClustersBestCandidate),
-                                    ptSpectrometer,
-                                    thetaSpectrometer,
-                                    etaSpectrometer);
-
-    fNtuFinalBestCandidates2 -> Fill(chi2HistoryForBestCandidate[0],
-                                    chi2HistoryForBestCandidate[1],
-                                    chi2HistoryForBestCandidate[2],
-                                    chi2HistoryForBestCandidate[3],
-                                    chi2HistoryForBestCandidate[4],
-                                    chi2HistoryForBestCandidate[5],
-                                    chi2HistoryForBestCandidate[6],
-                                    chi2HistoryForBestCandidate[7],
-                                    chi2HistoryForBestCandidate[8],
-                                    chi2HistoryForBestCandidate[9],
-                                    chi2HistoryForBestCandidate[10],
-                                    chi2HistoryForBestCandidate[11],
-                                    chi2HistoryForBestCandidate[12],
-                                    chi2HistoryForBestCandidate[13],
-                                    chi2HistoryForBestCandidate[14]);
-
+    newTrack -> SetNWrongClustersMC(newTrack->GetNMFTClusters() - nGoodClustersBestCandidate);
+    new ((*fMuonForwardTracks)[fMuonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack);
   }
+
+  //  fOutputQAFile->cd();
+  
+  Float_t finalBestCandidatesInfo[] = {Double_t(fRun),
+                                      Double_t(fEv),
+                                      Double_t(fCountRealTracksAnalyzedOfEvent),
+                                      Double_t(nFinalTracks),
+                                      Double_t(nClustersMC),
+                                      Double_t(nGoodClustersBestCandidate),
+                                      ptSpectrometer,
+                                      thetaSpectrometer,
+                                      etaSpectrometer,
+                                      chi2HistoryForBestCandidate[0],
+                                      chi2HistoryForBestCandidate[1],
+                                      chi2HistoryForBestCandidate[2],
+                                      chi2HistoryForBestCandidate[3],
+                                      chi2HistoryForBestCandidate[4],
+                                      chi2HistoryForBestCandidate[5],
+                                      chi2HistoryForBestCandidate[6],
+                                      chi2HistoryForBestCandidate[7],
+                                      chi2HistoryForBestCandidate[8],
+                                      nClustersPerPlane[0],
+                                      nClustersPerPlane[1],
+                                      nClustersPerPlane[2],
+                                      nClustersPerPlane[3],
+                                      nClustersPerPlane[4],
+                                      nClustersPerPlane[5],
+                                      nClustersPerPlane[6],
+                                      nClustersPerPlane[7],
+                                      nClustersPerPlane[8]};
+  
+  fNtuFinalBestCandidates -> Fill(finalBestCandidatesInfo);
   
   if (fDrawOption && bestCandidateExists) {
-    fTxtTrackGoodClusters = new TLatex(0.20, 0.59, Form("N_{GoodClusters} = %d", nGoodClustersBestCandidate));
+    fTxtTrackGoodClusters = new TLatex(0.20, 0.51, Form("N_{GoodClusters} = %d", nGoodClustersBestCandidate));
     DrawPlanes();
   }
 
   if (fIsCurrentMuonTrackable) {
+    //    fOutputQAFile->cd();
     if (nGoodClustersBestCandidate==5) fHistPtMuonTrackWithGoodMatch -> Fill(ptSpectrometer);
     else                               fHistPtMuonTrackWithBadMatch  -> Fill(ptSpectrometer);
   }
@@ -613,6 +773,7 @@ Int_t AliMuonForwardTrackFinder::LoadNextTrack() {
   // -------------------------------------------------------------------------------------------
 
   fCandidateTracks->Clear();
+  fFinalBestCandidate = NULL;
   
   fCountRealTracksAnalyzedOfEvent++;
 
@@ -624,7 +785,7 @@ Int_t AliMuonForwardTrackFinder::LoadNextTrack() {
 
 void AliMuonForwardTrackFinder::FindClusterInPlane(Int_t planeId) { 
   
-  printf(">>>> executing AliMuonForwardTrackFinder::FindClusterInPlane(%d)\n", 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
 
@@ -650,13 +811,13 @@ void AliMuonForwardTrackFinder::FindClusterInPlane(Int_t planeId) {
     currentParamForResearchFront = currentParamFront;
     currentParamForResearchBack  = currentParamBack;
     AliMUONTrackExtrap::AddMCSEffect(&currentParamFront,           (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
-                                                                   fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/radLengthSi,-1.);
+                                                                   fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
     AliMUONTrackExtrap::AddMCSEffect(&currentParamForResearchFront,(fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
-                                                                   fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/radLengthSi,-1.);
+                                                                   fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
     AliMUONTrackExtrap::AddMCSEffect(&currentParamBack,            (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
-                                                                   fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/radLengthSi,-1.);
+                                                                   fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
     AliMUONTrackExtrap::AddMCSEffect(&currentParamForResearchBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
-                                                                   fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/radLengthSi,-1.);
+                                                                   fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
   }
   // for all planes: extrapolation to the Z of the plane
   AliMUONTrackExtrap::ExtrapToZCov(&currentParamFront,            -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());   
@@ -681,7 +842,10 @@ void AliMuonForwardTrackFinder::FindClusterInPlane(Int_t planeId) {
   if (planeId==fNPlanesMFT-1 && 0.5*(researchRadiusFront+researchRadiusBack)<fMinResearchRadiusAtLastPlane) {
     corrFact = fMinResearchRadiusAtLastPlane/(0.5*(researchRadiusFront+researchRadiusBack));
   }
-  if (fIsCurrentMuonTrackable) fHistResearchRadius[planeId] -> Fill(corrFact*0.5*(researchRadiusFront+researchRadiusBack));
+  if (fIsCurrentMuonTrackable) {
+    //    fOutputQAFile->cd();
+    fHistResearchRadius[planeId] -> Fill(0.5*(researchRadiusFront+researchRadiusBack));
+  }
 
   Double_t position_X_Front = currentParamForResearchFront.GetNonBendingCoor();
   Double_t position_Y_Front = currentParamForResearchFront.GetBendingCoor();
@@ -697,7 +861,7 @@ void AliMuonForwardTrackFinder::FindClusterInPlane(Int_t planeId) {
   // Analyizing the clusters: FRONT ACTIVE ELEMENTS
   
   Int_t nClustersFront = fMFTClusterArrayFront[planeId]->GetEntries();
-  printf("There are %3d clusters in plane %02d FRONT\n", nClustersFront, planeId);
+  AliDebug(2, Form("There are %3d clusters in plane %02d FRONT\n", nClustersFront, planeId));
   
   for (Int_t iCluster=0; iCluster<nClustersFront; iCluster++) {
 
@@ -719,34 +883,38 @@ void AliMuonForwardTrackFinder::FindClusterInPlane(Int_t planeId) {
     }
 
     if (fIsCurrentMuonTrackable) {
+      //      fOutputQAFile->cd();
       if (IsCorrectMatch(cluster)) fHistChi2Cluster_GoodCluster[planeId]->Fill(chi2/2.);     //  chi2/ndf
       else                         fHistChi2Cluster_BadCluster[planeId] ->Fill(chi2/2.);     //  chi2/ndf
     }
 
     if (isGoodChi2) {
-      printf("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut);
+      AliDebug(3, Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
       AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
       newTrack->AddTrackParamAtMFTCluster(currentParamFront, *cluster);    // creating new track param and attaching the cluster
+      AliDebug(1, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)", 
+                      planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters()));
       newTrack->SetPlaneExists(planeId);
-      printf("current muon is trackable: %d\n", fIsCurrentMuonTrackable);
+      AliDebug(2, Form("current muon is trackable: %d\n", fIsCurrentMuonTrackable));
       if (fIsCurrentMuonTrackable) {
        Double_t newGlobalChi2 = ((AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First())->GetTrackChi2();
-       printf("new chi2 = %f (= %f)\n", newGlobalChi2, newTrack->GetMFTCluster(0)->GetTrackChi2());
+       AliDebug(2, Form("new chi2 = %f (= %f)\n", newGlobalChi2, newTrack->GetMFTCluster(0)->GetTrackChi2()));
        Int_t nClustersGlobalTrack = newTrack->GetNMUONClusters() + newTrack->GetNMFTClusters();        // Muon Spectrometer clusters + clusters in the Vertex Telescope
        Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
-       if (IsCorrectMatch(cluster)) fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[planeId]->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
-       else                         fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[planeId] ->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
+       //      fOutputQAFile->cd();
+       if (IsCorrectMatch(cluster)) fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[planeId]->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
+       else                         fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[planeId] ->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
       }
-      fGrMFTPlane[planeId][kClustersGoodChi2] -> SetPoint(fGrMFTPlane[planeId][kClustersGoodChi2]->GetN(), cluster->GetX(), cluster->GetY());
+      fGrMFTPlane[kClustersGoodChi2][planeId] -> SetPoint(fGrMFTPlane[kClustersGoodChi2][planeId]->GetN(), cluster->GetX(), cluster->GetY());
     }
-    else printf("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut);
+    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();
-  printf("There are %3d clusters in plane %02d BACK\n", nClustersBack, planeId);
+  AliDebug(2, Form("There are %3d clusters in plane %02d BACK\n", nClustersBack, planeId));
   
   for (Int_t iCluster=0; iCluster<nClustersBack; iCluster++) {
 
@@ -768,27 +936,31 @@ void AliMuonForwardTrackFinder::FindClusterInPlane(Int_t planeId) {
     }
 
     if (fIsCurrentMuonTrackable) {
+      //      fOutputQAFile->cd();
       if (IsCorrectMatch(cluster)) fHistChi2Cluster_GoodCluster[planeId]->Fill(chi2/2.);     //  chi2/ndf
       else                         fHistChi2Cluster_BadCluster[planeId] ->Fill(chi2/2.);     //  chi2/ndf
     }
 
     if (isGoodChi2) {
-      printf("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut);
+      AliDebug(3,Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
       AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
       newTrack->AddTrackParamAtMFTCluster(currentParamBack, *cluster);    // creating new track param and attaching the cluster
+      AliDebug(1, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)", 
+                      planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters()));
       newTrack->SetPlaneExists(planeId);
-      printf("current muon is trackable: %d\n", fIsCurrentMuonTrackable);
+      AliDebug(2, Form("current muon is trackable: %d\n", fIsCurrentMuonTrackable));
       if (fIsCurrentMuonTrackable) {
        Double_t newGlobalChi2 = ((AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First())->GetTrackChi2();
-       printf("new chi2 = %f (= %f)\n", newGlobalChi2, newTrack->GetMFTCluster(0)->GetTrackChi2());
+       AliDebug(2, Form("new chi2 = %f (= %f)\n", newGlobalChi2, newTrack->GetMFTCluster(0)->GetTrackChi2()));
        Int_t nClustersGlobalTrack = newTrack->GetNMUONClusters() + newTrack->GetNMFTClusters();        // Muon Spectrometer clusters + clusters in the Vertex Telescope
        Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
-       if (IsCorrectMatch(cluster)) fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[planeId]->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
-       else                         fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[planeId] ->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
+       //      fOutputQAFile->cd();
+       if (IsCorrectMatch(cluster)) fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[planeId]->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
+       else                         fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[planeId] ->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
       }
-      fGrMFTPlane[planeId][kClustersGoodChi2] -> SetPoint(fGrMFTPlane[planeId][kClustersGoodChi2]->GetN(), cluster->GetX(), cluster->GetY());
+      fGrMFTPlane[kClustersGoodChi2][planeId] -> SetPoint(fGrMFTPlane[kClustersGoodChi2][planeId]->GetN(), cluster->GetX(), cluster->GetY());
     }
-    else printf("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut);
+    else AliDebug(3,Form("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
 
   }
 
@@ -796,6 +968,7 @@ void AliMuonForwardTrackFinder::FindClusterInPlane(Int_t planeId) {
 
   if (planeId == fNPlanesMFT-1) {
     if (fIsCurrentMuonTrackable && fDistanceFromGoodClusterAndTrackAtLastPlane>0.) {
+      //      fOutputQAFile->cd();
       fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> Fill(TMath::Abs(fDistanceFromBestClusterAndTrackAtLastPlane-
                                                                                                       fDistanceFromGoodClusterAndTrackAtLastPlane));
       fHistDistanceGoodClusterFromTrackAtLastPlane -> Fill(fDistanceFromGoodClusterAndTrackAtLastPlane);
@@ -808,7 +981,7 @@ void AliMuonForwardTrackFinder::FindClusterInPlane(Int_t planeId) {
 
 void AliMuonForwardTrackFinder::AttachGoodClusterInPlane(Int_t planeId) { 
   
-  printf(">>>> executing AliMuonForwardTrackFinder::AttachGoodClusterInPlane(%d)\n", planeId);
+  AliDebug(1, Form(">>>> executing AliMuonForwardTrackFinder::AttachGoodClusterInPlane(%d)\n", planeId));
 
   AliMUONTrackParam currentParamFront, currentParamBack;
 
@@ -819,13 +992,13 @@ void AliMuonForwardTrackFinder::AttachGoodClusterInPlane(Int_t planeId) {
     AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamBack,  0.); 
   }
   else {          // MFT planes others than the last one: mult. scattering correction because of the upstream MFT planes is performed
-    printf("fCurrentTrack = %p\n", fCurrentTrack);
+    AliDebug(2, Form("fCurrentTrack = %p\n", fCurrentTrack));
     currentParamFront = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
     currentParamBack  = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
     AliMUONTrackExtrap::AddMCSEffect(&currentParamFront, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
-                                                         fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/radLengthSi,-1.);
+                                                         fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
     AliMUONTrackExtrap::AddMCSEffect(&currentParamBack,  (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
-                                                         fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/radLengthSi,-1.);
+                                                         fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
   }
   // for all planes: linear extrapolation to the Z of the plane
   AliMUONTrackExtrap::ExtrapToZCov(&currentParamFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());   
@@ -837,10 +1010,10 @@ void AliMuonForwardTrackFinder::AttachGoodClusterInPlane(Int_t planeId) {
 
   Int_t nClustersFront = fMFTClusterArrayFront[planeId]->GetEntries();
   
-  printf("nClustersFront = %d\n", nClustersFront);
+  AliDebug(1, Form("nClustersFront = %d\n", nClustersFront));
   for (Int_t iCluster=0; iCluster<nClustersFront; iCluster++) {
     AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayFront[planeId]->UncheckedAt(iCluster);
-    printf("checking cluster %02d of %02d: cluter=%p, fCurrentTrack=%p\n", iCluster, nClustersFront, cluster, fCurrentTrack);
+    AliDebug(2, Form("checking cluster %02d of %02d: cluter=%p, fCurrentTrack=%p\n", iCluster, nClustersFront, cluster, fCurrentTrack));
     if (IsCorrectMatch(cluster)) {
       fCurrentTrack->AddTrackParamAtMFTCluster(currentParamFront, *cluster);  // creating new track param and attaching the cluster
       fCurrentTrack->SetPlaneExists(planeId);
@@ -855,10 +1028,10 @@ void AliMuonForwardTrackFinder::AttachGoodClusterInPlane(Int_t planeId) {
 
   Int_t nClustersBack = fMFTClusterArrayBack[planeId]->GetEntries();
   
-  printf("nClustersBack = %d\n", nClustersBack);
+  AliDebug(1, Form("nClustersBack = %d\n", nClustersBack));
   for (Int_t iCluster=0; iCluster<nClustersBack; iCluster++) {
     AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayBack[planeId]->UncheckedAt(iCluster);
-    printf("checking cluster %02d of %02d: cluter=%p, fCurrentTrack=%p\n", iCluster, nClustersBack, cluster, fCurrentTrack);
+    AliDebug(2,Form("checking cluster %02d of %02d: cluter=%p, fCurrentTrack=%p\n", iCluster, nClustersBack, cluster, fCurrentTrack));
     if (IsCorrectMatch(cluster)) {
       fCurrentTrack->AddTrackParamAtMFTCluster(currentParamBack, *cluster);  // creating new track param and attaching the cluster
       fCurrentTrack->SetPlaneExists(planeId);
@@ -894,30 +1067,25 @@ void AliMuonForwardTrackFinder::CheckCurrentMuonTrackable() {
 
 //==========================================================================================================================================
 
-void AliMuonForwardTrackFinder::FillPlanesWithTrackHistory(AliMuonForwardTrack *track) { 
+void AliMuonForwardTrackFinder::FillPlanesWithTrackHistory() { 
   
-  // recover track parameters on each planes and look for the corresponding clusters
+  // Fill planes with the clusters
 
+  Int_t cluster = 0;
+  AliDebug(2, Form("fFinalBestCandidate->GetNMFTClusters() = %d\n",  fFinalBestCandidate->GetNMFTClusters()));
   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
-    
-    AliMFTCluster *trackCluster = (AliMFTCluster *) track->GetMFTCluster(iPlane);
-
-    fGrMFTPlane[iPlane][kClusterOfTrack] -> SetPoint(fGrMFTPlane[iPlane][kClusterOfTrack]->GetN(), trackCluster->GetX(), trackCluster->GetY());
-    
+    if (fFinalBestCandidate->PlaneExists(iPlane)) {
+      AliMFTCluster *trackCluster = fFinalBestCandidate->GetMFTCluster(cluster++);
+      fGrMFTPlane[kClusterOfTrack][iPlane] -> SetPoint(fGrMFTPlane[kClusterOfTrack][iPlane]->GetN(), trackCluster->GetX(), trackCluster->GetY());
+    }
     Int_t nClusters = fMFTClusterArray[iPlane]->GetEntriesFast();
-    
     for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
-      
-      AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->UncheckedAt(iCluster); 
-
-      fGrMFTPlane[iPlane][kAllClusters] -> SetPoint(fGrMFTPlane[iPlane][kAllClusters]->GetN(), cluster->GetX(), cluster->GetY());
-      
-      if (IsCorrectMatch(cluster)) {
-       fGrMFTPlane[iPlane][kClusterCorrectMC] -> SetPoint(fGrMFTPlane[iPlane][kClusterCorrectMC]->GetN(), cluster->GetX(), cluster->GetY());
+      AliMFTCluster *myCluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->UncheckedAt(iCluster); 
+      fGrMFTPlane[kAllClusters][iPlane] -> SetPoint(fGrMFTPlane[kAllClusters][iPlane]->GetN(), myCluster->GetX(), myCluster->GetY());
+      if (IsCorrectMatch(myCluster)) {
+       fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetPoint(fGrMFTPlane[kClusterCorrectMC][iPlane]->GetN(), myCluster->GetX(), myCluster->GetY());
       }
-
     }
-
   }
 
 }
@@ -937,7 +1105,7 @@ Bool_t AliMuonForwardTrackFinder::IsCorrectMatch(AliMFTCluster *cluster) {
     }
   }
 
-  printf("returning %d\n", result);
+  AliDebug(2,Form("returning %d\n", result));
 
   return result;
 
@@ -954,13 +1122,13 @@ Double_t AliMuonForwardTrackFinder::TryOneCluster(const AliMUONTrackParam &track
   // 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();
-  printf("dX = %f, dY = %f\n", dX, dY);
+  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();
-  printf("dX2 = %f, dY2 = %f\n", sigmaX2, sigmaY2);
+  AliDebug(3, Form("dX2 = %f, dY2 = %f\n", sigmaX2, sigmaY2));
   Double_t covXY   = kParamCov(0,2);
   Double_t det     = sigmaX2 * sigmaY2 - covXY * covXY;
   
@@ -1004,9 +1172,9 @@ Int_t AliMuonForwardTrackFinder::GetNDF(Int_t nClusters) {
 //============================================================================================================================================
 
 void AliMuonForwardTrackFinder::BookHistos() {
-  
+
   const Int_t nMaxNewTracks[]  = {150,     200,   250, 600, 1000};
-  const Double_t radiusPlane[] = {0.001, 0.010, 0.100, 5.0,  5.0};
+  const Double_t radiusPlane[] = {0.010, 0.010, 0.050, 0.5,  1.5};
 
   fHistPtSpectrometer = new TH1D("hPtSpectrometer", "p_{T} as given by the Muon Spectrometer", 200, 0, 20.); 
 
@@ -1042,14 +1210,14 @@ void AliMuonForwardTrackFinder::BookHistos() {
                                                   Form("#chi^{2}_{clust} for Bad clusters in MFT plane %02d", iPlane),
                                                   100, 0., 15.);
 
-    fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] = new TH1D(Form("fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons_pl%02d", iPlane),
-                                                                          Form("#chi^{2}/ndf at plane %d for GOOD candidates of trackable muons",iPlane),
-                                                                          100, 0., 15.);
-
-    fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] = new TH1D(Form("fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons_pl%02d", iPlane),
-                                                                          Form("#chi^{2}/ndf at plane %d for BAD candidates of trackable muons",iPlane),
-                                                                          100, 0., 15.);
-
+    fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] = new TH1D(Form("fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons_pl%02d", iPlane),
+                                                                                Form("#chi^{2}/ndf at plane %d for GOOD candidates of trackable muons",iPlane),
+                                                                                100, 0., 15.);
+    
+    fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] = new TH1D(Form("fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons_pl%02d", iPlane),
+                                                                               Form("#chi^{2}/ndf at plane %d for BAD candidates of trackable muons",iPlane),
+                                                                               100, 0., 15.);
+    
   }
   
   //------------------------------------------
@@ -1071,18 +1239,14 @@ void AliMuonForwardTrackFinder::BookHistos() {
     fHistChi2Cluster_GoodCluster[iPlane]        -> Sumw2();
     fHistChi2Cluster_BadCluster[iPlane]         -> Sumw2();
 
-    fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> Sumw2();
-    fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane]  -> Sumw2();
+    fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> Sumw2();
+    fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane]  -> Sumw2();
     
   }
 
-  fNtuFinalCandidates1     = new TNtuple("ntuFinalCandidates1",     "Final Candidates (ALL)", "run:event:muonTrack:nFinalCandidates:nClustersMC:nGoodClusters:ptSpectrometer:thetaSpectrometer:etaSpectrometer");
-
-  fNtuFinalBestCandidates1 = new TNtuple("ntuFinalBestCandidates1", "Final Best Candidates",  "run:event:muonTrack:nFinalCandidates:nClustersMC:nGoodClusters:ptSpectrometer:thetaSpectrometer:etaSpectrometer");
+  fNtuFinalCandidates     = new TNtuple("ntuFinalCandidates",     "Final Candidates (ALL)", "run:event:muonTrack:nFinalCandidates:nClustersMC:nGoodClusters:ptSpectrometer:thetaSpectrometer:etaSpectrometer:chi2AtPlane0:chi2AtPlane1:chi2AtPlane2:chi2AtPlane3:chi2AtPlane4:chi2AtPlane5:chi2AtPlane6:chi2AtPlane7:chi2AtPlane8");
 
-  fNtuFinalCandidates2     = new TNtuple("ntuFinalCandidates2",     "Final Candidates (ALL)", "chi2AtPlane0:chi2AtPlane1:chi2AtPlane2:chi2AtPlane3:chi2AtPlane4:chi2AtPlane5:chi2AtPlane6:chi2AtPlane7:chi2AtPlane8:chi2AtPlane9:chi2AtPlane10:chi2AtPlane11:chi2AtPlane12:chi2AtPlane13:chi2AtPlane14");
-
-  fNtuFinalBestCandidates2 = new TNtuple("ntuFinalBestCandidates2", "Final Best Candidates",  "chi2AtPlane0:chi2AtPlane1:chi2AtPlane2:chi2AtPlane3:chi2AtPlane4:chi2AtPlane5:chi2AtPlane6:chi2AtPlane7:chi2AtPlane8:chi2AtPlane9:chi2AtPlane10:chi2AtPlane11:chi2AtPlane12:chi2AtPlane13:chi2AtPlane14");
+  fNtuFinalBestCandidates = new TNtuple("ntuFinalBestCandidates", "Final Best Candidates",  "run:event:muonTrack:nFinalCandidates:nClustersMC:nGoodClusters:ptSpectrometer:thetaSpectrometer:etaSpectrometer:chi2AtPlane0:chi2AtPlane1:chi2AtPlane2:chi2AtPlane3:chi2AtPlane4:chi2AtPlane5:chi2AtPlane6:chi2AtPlane7:chi2AtPlane8:nClustersAtPlane0:nClustersAtPlane1:nClustersAtPlane2:nClustersAtPlane3:nClustersAtPlane4:nClustersAtPlane5:nClustersAtPlane6:nClustersAtPlane7:nClustersAtPlane8");
 
 }
 
@@ -1108,8 +1272,8 @@ void AliMuonForwardTrackFinder::SetTitleHistos() {
     fHistChi2Cluster_GoodCluster[iPlane]         -> SetXTitle("#chi^{2}/ndf");
     fHistChi2Cluster_BadCluster[iPlane]          -> SetXTitle("#chi^{2}/ndf");
 
-    fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> SetXTitle("#chi^{2}/ndf");
-    fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane]  -> SetXTitle("#chi^{2}/ndf");
+    fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> SetXTitle("#chi^{2}/ndf");
+    fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane]  -> SetXTitle("#chi^{2}/ndf");
     
   }
 
@@ -1120,41 +1284,41 @@ void AliMuonForwardTrackFinder::SetTitleHistos() {
 void AliMuonForwardTrackFinder::BookPlanes() {
 
   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
-    fGrMFTPlane[iPlane][kAllClusters] = new TGraph();
-    fGrMFTPlane[iPlane][kAllClusters] -> SetName(Form("fGrMFTPlane_%02d_AllClusters",iPlane));
-    fGrMFTPlane[iPlane][kAllClusters] -> SetMarkerStyle(20);
-    //    fGrMFTPlane[iPlane][kAllClusters] -> SetMarkerSize(0.5);
-    //    fGrMFTPlane[iPlane][kAllClusters] -> SetMarkerSize(0.3);
-    fGrMFTPlane[iPlane][kAllClusters] -> SetMarkerSize(0.2);
+    fGrMFTPlane[kAllClusters][iPlane] = new TGraph();
+    fGrMFTPlane[kAllClusters][iPlane] -> SetName(Form("fGrMFTPlane_%02d_AllClusters",iPlane));
+    fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerStyle(20);
+    //    fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerSize(0.5);
+    //    fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerSize(0.3);
+    fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerSize(0.2);
   }
 
   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
-    fGrMFTPlane[iPlane][kClustersGoodChi2] = new TGraph();
-    fGrMFTPlane[iPlane][kClustersGoodChi2] -> SetName(Form("fGrMFTPlane_%02d_ClustersGoodChi2",iPlane));
-    fGrMFTPlane[iPlane][kClustersGoodChi2] -> SetMarkerStyle(20);
-    //    fGrMFTPlane[iPlane][kClustersGoodChi2] -> SetMarkerSize(0.8);
-    //    fGrMFTPlane[iPlane][kClustersGoodChi2] -> SetMarkerSize(0.4);
-    fGrMFTPlane[iPlane][kClustersGoodChi2] -> SetMarkerSize(0.3);
-    fGrMFTPlane[iPlane][kClustersGoodChi2] -> SetMarkerColor(kBlue);
+    fGrMFTPlane[kClustersGoodChi2][iPlane] = new TGraph();
+    fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetName(Form("fGrMFTPlane_%02d_ClustersGoodChi2",iPlane));
+    fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerStyle(20);
+    //    fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerSize(0.8);
+    //    fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerSize(0.4);
+    fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerSize(0.3);
+    fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerColor(kBlue);
   }
 
   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
-    fGrMFTPlane[iPlane][kClusterOfTrack] = new TGraph();
-    fGrMFTPlane[iPlane][kClusterOfTrack] -> SetName(Form("fGrMFTPlane_%02d_ClustersOfTrack",iPlane));
-    fGrMFTPlane[iPlane][kClusterOfTrack] -> SetMarkerStyle(25);
-    //    fGrMFTPlane[iPlane][kClusterOfTrack] -> SetMarkerSize(1.2);
-    fGrMFTPlane[iPlane][kClusterOfTrack] -> SetMarkerSize(0.9);
-    fGrMFTPlane[iPlane][kClusterOfTrack] -> SetMarkerColor(kRed);
-    fGrMFTPlane[iPlane][kClusterOfTrack] -> SetTitle(Form("Plane %d (%3.1f cm)", iPlane, fZPlane[iPlane]));
+    fGrMFTPlane[kClusterOfTrack][iPlane] = new TGraph();
+    fGrMFTPlane[kClusterOfTrack][iPlane] -> SetName(Form("fGrMFTPlane_%02d_ClustersOfTrack",iPlane));
+    fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerStyle(25);
+    //    fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerSize(1.2);
+    fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerSize(0.9);
+    fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerColor(kRed);
+    fGrMFTPlane[kClusterOfTrack][iPlane] -> SetTitle(Form("Plane %d (%3.1f cm)", iPlane, fZPlane[iPlane]));
   }
 
   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
-    fGrMFTPlane[iPlane][kClusterCorrectMC] = new TGraph();
-    fGrMFTPlane[iPlane][kClusterCorrectMC] -> SetName(Form("fGrMFTPlane_%02d_ClustersCorrectMC",iPlane));
-    fGrMFTPlane[iPlane][kClusterCorrectMC] -> SetMarkerStyle(20);
-    //    fGrMFTPlane[iPlane][kClusterCorrectMC] -> SetMarkerSize(0.8);
-    fGrMFTPlane[iPlane][kClusterCorrectMC] -> SetMarkerSize(0.5);
-    fGrMFTPlane[iPlane][kClusterCorrectMC] -> SetMarkerColor(kGreen);
+    fGrMFTPlane[kClusterCorrectMC][iPlane] = new TGraph();
+    fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetName(Form("fGrMFTPlane_%02d_ClustersCorrectMC",iPlane));
+    fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerStyle(20);
+    //    fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerSize(0.8);
+    fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerSize(0.5);
+    fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerColor(kGreen);
   }
 
   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
@@ -1162,7 +1326,7 @@ void AliMuonForwardTrackFinder::BookPlanes() {
     fCircleInt[iPlane] = new TEllipse(0., 0., fRPlaneMin[iPlane], fRPlaneMin[iPlane]);
   }
   
-  fTxtDummy = new TLatex(0.10, 0.67, "Best Candidate:");
+  fTxtDummy = new TLatex(0.10, 0.59, "Best Candidate:");
 
   //---------------------------------------------------
 
@@ -1201,8 +1365,8 @@ void AliMuonForwardTrackFinder::ResetPlanes() {
 
   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
     for (Int_t iGr=0; iGr<4; iGr++) {
-      Int_t nOldClusters = fGrMFTPlane[iPlane][iGr]->GetN();
-      for (Int_t iPoint=nOldClusters-1; iPoint>=0; iPoint--) fGrMFTPlane[iPlane][iGr]->RemovePoint(iPoint);
+      Int_t nOldClusters = fGrMFTPlane[iGr][iPlane]->GetN();
+      for (Int_t iPoint=nOldClusters-1; iPoint>=0; iPoint--) fGrMFTPlane[iGr][iPlane]->RemovePoint(iPoint);
     }
   }
 
@@ -1212,38 +1376,66 @@ void AliMuonForwardTrackFinder::ResetPlanes() {
 
 void AliMuonForwardTrackFinder::PrintParticleHistory() {
 
+  AliDebug(1, "Entering");
+
   TString history = "";
   
-  TParticle *part = fStack->Particle(fLabelMC);
-  
-  if (part->GetFirstMother() != -1) {
-    TParticle *partMother = fStack->Particle(part->GetFirstMother());
-    if (partMother->GetFirstMother() != -1) history += "...  #rightarrow ";
+  TParticle *part = 0;
+  if (fLabelMC>=0) part = fStack->Particle(fLabelMC);
+
+  AliDebug(1, Form("fStack->Particle(fLabelMC) = %p", part));
+
+  if (part) {
+    if (part->GetFirstMother() != -1) {
+      TParticle *partMother = fStack->Particle(part->GetFirstMother());
+      AliDebug(1, Form("fStack->Particle(part->GetFirstMother() = %p", partMother));
+      if (partMother) {
+       Char_t newName[100];
+       if (partMother->GetFirstMother() != -1) history += "...  #rightarrow ";
+       PDGNameConverter(partMother->GetName(), newName);
+       history += Form("%s #rightarrow ", newName);
+      }
+    }
     Char_t newName[100];
-    PDGNameConverter(partMother->GetName(), newName);
-    history += Form("%s #rightarrow ", newName);
+    PDGNameConverter(part->GetName(), newName);
+    history += Form("%s  at  z = %5.1f cm", newName, part->Vz());
+    //  printf("%s", history.Data());
   }
-  Char_t newName[100];
-  PDGNameConverter(part->GetName(), newName);
-  history += Form("%s  at  z = %5.1f cm", newName, part->Vz());
-  
-  //  printf("%s", history.Data());
-  
+  else history += "NO AVAILABLE HISTORY";
+
   fTxtMuonHistory = new TLatex(0.10, 0.86, history.Data());
+
+  // Filling particle history in the fFinalBestCandidate
+
+  if (part) {
+    for (Int_t iParent=0; iParent<AliMuonForwardTrack::fgkNParentsMax; iParent++) {
+      if (part->GetFirstMother() == -1) break;
+      if (!(fStack->Particle(part->GetFirstMother()))) break;
+      AliDebug(1, Form("fStack->Particle(part->GetFirstMother() = %p", fStack->Particle(part->GetFirstMother())));
+      fFinalBestCandidate->SetParentMCLabel(iParent, part->GetFirstMother());
+      fFinalBestCandidate->SetParentPDGCode(iParent, fStack->Particle(part->GetFirstMother())->GetPdgCode());
+      part = fStack->Particle(part->GetFirstMother());
+    }
+  }
   
 }
 
 //===========================================================================================================================================
 
-Bool_t AliMuonForwardTrackFinder::IsMother(const Char_t *nameMother) {
+Bool_t AliMuonForwardTrackFinder::IsMother(Char_t *nameMother) {
   
   Bool_t result = kFALSE;
   
-  TParticle *part = fStack->Particle(fLabelMC);
+  TParticle *part = 0;
+  if (fLabelMC>=0) part = fStack->Particle(fLabelMC);
   
-  if (part->GetFirstMother() != -1) {
-    TParticle *partMother = fStack->Particle(part->GetFirstMother());
-    if (!strcmp(partMother->GetName(), nameMother)) result=kTRUE;
+  if (part) {
+    if (part->GetFirstMother() != -1) {
+      TParticle *partMother = fStack->Particle(part->GetFirstMother());
+      if (partMother) {
+       if (!strcmp(partMother->GetName(), nameMother)) result=kTRUE;
+      }
+    }
   }
 
   return result;
@@ -1263,19 +1455,19 @@ void AliMuonForwardTrackFinder::DrawPlanes() {
     
     fCanvas->cd(fNPlanesMFT-iPlane+1);
     
-    fGrMFTPlane[iPlane][kClusterOfTrack] -> GetXaxis() -> SetLimits(-1.1*fRPlaneMax[fNPlanesMFT-1], +1.1*fRPlaneMax[fNPlanesMFT-1]);
-    fGrMFTPlane[iPlane][kClusterOfTrack] -> GetYaxis() -> SetRangeUser(-1.1*fRPlaneMax[fNPlanesMFT-1], +1.1*fRPlaneMax[fNPlanesMFT-1]);
-    fGrMFTPlane[iPlane][kClusterOfTrack] -> GetXaxis() -> SetTitle("X  [cm]");
-    fGrMFTPlane[iPlane][kClusterOfTrack] -> GetYaxis() -> SetTitle("Y  [cm]");
-    fGrMFTPlane[iPlane][kClusterOfTrack] -> Draw("ap");
+    fGrMFTPlane[kClusterOfTrack][iPlane] -> GetXaxis() -> SetLimits(-1.1*fRPlaneMax[fNPlanesMFT-1], +1.1*fRPlaneMax[fNPlanesMFT-1]);
+    fGrMFTPlane[kClusterOfTrack][iPlane] -> GetYaxis() -> SetRangeUser(-1.1*fRPlaneMax[fNPlanesMFT-1], +1.1*fRPlaneMax[fNPlanesMFT-1]);
+    fGrMFTPlane[kClusterOfTrack][iPlane] -> GetXaxis() -> SetTitle("X  [cm]");
+    fGrMFTPlane[kClusterOfTrack][iPlane] -> GetYaxis() -> SetTitle("Y  [cm]");
+    fGrMFTPlane[kClusterOfTrack][iPlane] -> Draw("ap");
 
     fCircleExt[iPlane] -> Draw("same");
     fCircleInt[iPlane] -> Draw("same");
     
-    if (fGrMFTPlane[iPlane][kAllClusters]->GetN())       fGrMFTPlane[iPlane][kAllClusters]      -> Draw("psame");
-    if (fGrMFTPlane[iPlane][kClustersGoodChi2]->GetN())  fGrMFTPlane[iPlane][kClustersGoodChi2] -> Draw("psame");
-    if (fGrMFTPlane[iPlane][kClusterOfTrack]->GetN())    fGrMFTPlane[iPlane][kClusterOfTrack]   -> Draw("psame");
-    if (fGrMFTPlane[iPlane][kClusterCorrectMC]->GetN())  fGrMFTPlane[iPlane][kClusterCorrectMC] -> Draw("psame");
+    if (fGrMFTPlane[kAllClusters][iPlane]->GetN())       fGrMFTPlane[kAllClusters][iPlane]      -> Draw("psame");
+    if (fGrMFTPlane[kClustersGoodChi2][iPlane]->GetN())  fGrMFTPlane[kClustersGoodChi2][iPlane] -> Draw("psame");
+    if (fGrMFTPlane[kClusterOfTrack][iPlane]->GetN())    fGrMFTPlane[kClusterOfTrack][iPlane]   -> Draw("psame");
+    if (fGrMFTPlane[kClusterCorrectMC][iPlane]->GetN())  fGrMFTPlane[kClusterCorrectMC][iPlane] -> Draw("psame");
 
     fTxtTrackChi2[iPlane] -> Draw("same");
 
@@ -1284,9 +1476,10 @@ void AliMuonForwardTrackFinder::DrawPlanes() {
   fCanvas -> cd(1);
   fTxtMuonHistory       -> Draw();
   fTxtDummy             -> Draw("same");
-  fTxtTrackGoodClusters -> Draw("same");
+  if (fMatchingMode==kRealMatching) fTxtTrackGoodClusters -> Draw("same");
   fTxtTrackFinalChi2    -> Draw("same");
-  fTxtFinalCandidates   -> Draw("same");
+  fTxtTrackMomentum     -> Draw("same");
+  if (fMatchingMode==kRealMatching) fTxtFinalCandidates   -> Draw("same");
 
   fMrkAllClust      -> Draw("same");
   fMrkClustGoodChi2 -> Draw("same");
@@ -1315,15 +1508,20 @@ void AliMuonForwardTrackFinder::DrawPlanes() {
 
 void AliMuonForwardTrackFinder::Terminate() {
   
-  printf("\n");
-  printf("---------------------------------------------------------------------------------------------------------------\n");
-  printf("%8d  tracks analyzed\n",                                                     fCountRealTracksAnalyzed);
-  printf("%8d  tracks with MC ref\n",                                                  fCountRealTracksWithRefMC);
-  printf("%8d  tracks with MC ref & trigger match\n",                                  fCountRealTracksWithRefMC_andTrigger);
-  printf("%8d  tracks analyzed with final candidates\n",                               fCountRealTracksAnalyzedWithFinalCandidates);
-//   printf("%8d  tracks with MC ref & trigger match & pt>%3.1f GeV/c\n",                 fCountRealTracksWithRefMC_andTrigger_andGoodPt, fLowPtCut);
-//   printf("%8d  tracks with MC ref & trigger match & pt>%3.1f GeV/c & correct R_abs\n", fCountRealTracksWithRefMC_andTrigger_andGoodPt_andGoodTheta, fLowPtCut);
-  printf("---------------------------------------------------------------------------------------------------------------\n");
+  AliInfo("");
+  AliInfo("---------------------------------------------------------------------------------------------------------------");
+  AliInfo(Form("%8d  tracks analyzed",                                                     fCountRealTracksAnalyzed));
+  AliInfo(Form("%8d  tracks with MC ref",                                                  fCountRealTracksWithRefMC));
+  AliInfo(Form("%8d  tracks with MC ref & trigger match",                                  fCountRealTracksWithRefMC_andTrigger));
+  if (fMatchingMode==kRealMatching) {
+    AliInfo(Form("%8d  tracks analyzed with final candidates",                             fCountRealTracksAnalyzedWithFinalCandidates));
+  }
+  else {
+    AliInfo(Form("%8d  tracks matched with their MC clusters",                             fCountRealTracksAnalyzedWithFinalCandidates));
+  }
+//   printf("%8d  tracks with MC ref & trigger match & pt>%3.1f GeV/c",                 fCountRealTracksWithRefMC_andTrigger_andGoodPt, fLowPtCut);
+//   printf("%8d  tracks with MC ref & trigger match & pt>%3.1f GeV/c & correct R_abs", fCountRealTracksWithRefMC_andTrigger_andGoodPt_andGoodTheta, fLowPtCut);
+  AliInfo("---------------------------------------------------------------------------------------------------------------");
 
   WriteOutputTree();
   WriteHistos();
@@ -1339,8 +1537,9 @@ void AliMuonForwardTrackFinder::FillOutputTree() {
   AliDebug(1, Form("Filling output tree %p with %p having %d entries whose 1st entry is %p", 
                   fOutputEventTree, fMuonForwardTracks, fMuonForwardTracks->GetEntries(), fMuonForwardTracks->At(0)));
   
+  //  fOutputTreeFile->cd();
   fOutputEventTree->Fill();
-  AliDebug(1, Form("\n\nFilled Tree: nEvents = %d!!!!\n\n", Int_t(fOutputEventTree->GetEntries())));
+  AliDebug(1, Form("\nFilled Tree: nEvents = %d!!!!\n", Int_t(fOutputEventTree->GetEntries())));
 
 }
 
@@ -1361,7 +1560,8 @@ void AliMuonForwardTrackFinder::WriteOutputTree() {
 
 void AliMuonForwardTrackFinder::WriteHistos() {
 
-  TFile *fileOut = new TFile(Form("%s/MuonGlobalTracking.QA.run%d.root", fOutDir.Data(), fRun), "recreate");
+  fOutputQAFile = new TFile(Form("MuonGlobalTracking.QA.run%d.root", fRun), "recreate");
+  fOutputQAFile -> cd();
 
   fHistPtSpectrometer              -> Write();
   fHistPtMuonTrackWithGoodMatch    -> Write();
@@ -1380,17 +1580,15 @@ void AliMuonForwardTrackFinder::WriteHistos() {
     fHistChi2Cluster_GoodCluster[iPlane]   -> Write();
     fHistChi2Cluster_BadCluster[iPlane]    -> Write();
     
-    fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> Write();
-    fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane]  -> Write();
+    fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> Write();
+    fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane]  -> Write();
 
   }
 
-  fNtuFinalCandidates1     -> Write();
-  fNtuFinalBestCandidates1 -> Write();
-  fNtuFinalCandidates2     -> Write();
-  fNtuFinalBestCandidates2 -> Write();
+  fNtuFinalCandidates     -> Write();
+  fNtuFinalBestCandidates -> Write();
 
-  fileOut -> Close();
+  fOutputQAFile -> Close();
 
 }
 
index 3cc512c..c26da03 100644 (file)
@@ -48,6 +48,7 @@
 #include "AliMFTCluster.h"
 #include "AliMFT.h"
 #include "AliMFTSegmentation.h"
+#include "AliMFTConstants.h"
 
 //====================================================================================================================================================
 //
@@ -66,7 +67,7 @@ public:
 
   AliMuonForwardTrackFinder();
   
-  virtual ~AliMuonForwardTrackFinder() {;}
+  virtual ~AliMuonForwardTrackFinder();
   
   enum {kAllClusters, kClustersGoodChi2, kClusterOfTrack, kClusterCorrectMC};
 
@@ -115,11 +116,11 @@ public:
   void SetNPlanesMFT(Int_t nPlanesMFT) { fNPlanesMFT = nPlanesMFT; }
   void SeparateFrontBackClusters();
   void SetNMaxMissingMFTClusters(Int_t nMaxMissingMFTClusters) { fNMaxMissingMFTClusters = nMaxMissingMFTClusters; }
-  void SetMandatoryPlane(Int_t iPlane) { if (0<=iPlane && iPlane<fMaxNPlanesMFT) fIsPlaneMandatory[iPlane] = kTRUE; }
+  void SetMandatoryPlane(Int_t iPlane) { if (0<=iPlane && iPlane<AliMFTConstants::fNMaxPlanes) fIsPlaneMandatory[iPlane] = kTRUE; }
 
   void FindClusterInPlane(Int_t planeId);
   void AttachGoodClusterInPlane(Int_t planeId);
-  void FillPlanesWithTrackHistory(AliMuonForwardTrack *track);
+  void FillPlanesWithTrackHistory();
   Double_t TryOneCluster(const AliMUONTrackParam &trackParam, AliMFTCluster *cluster);
   void BookHistos();
   void SetTitleHistos();
@@ -129,7 +130,7 @@ public:
 
   Bool_t IsCorrectMatch(AliMFTCluster *cluster);
   void CheckCurrentMuonTrackable();
-  Bool_t IsMother(const Char_t *nameMother);
+  Bool_t IsMother(Char_t *nameMother);
 
   void SetMatchingMode(Int_t matchingMode) { fMatchingMode = matchingMode; }
   void SetMinResearchRadiusAtLastPlane(Double_t minResearchRadius) { fMinResearchRadiusAtLastPlane = minResearchRadius; }
@@ -140,6 +141,8 @@ public:
   Bool_t InitGRP();
   Bool_t SetRunNumber();
 
+  void SetMaxNTracksToBeAnalyzed(Int_t nTracks) { fMaxNTracksToBeAnalyzed = nTracks; }
+
 private:
 
   AliMuonForwardTrackFinder(const AliMuonForwardTrackFinder& obj);
@@ -147,9 +150,9 @@ private:
 
 protected:
 
-  static const Int_t fMaxNPlanesMFT = 20;
-  static const Double_t radLengthSi = 9.37;  // expressed in cm
-  
+  static const Int_t fNMaxPlanes = AliMFTConstants::fNMaxPlanes;        // max number of MFT planes
+  static const Double_t fRadLengthSi;
+
   Int_t fRun; 
   Int_t fNEventsToAnalyze;           // events to analyze
   Double_t fSigmaClusterCut;         // to select the clusters in the MFT planes which are compatible with the extrapolated muon track
@@ -165,46 +168,49 @@ protected:
   Double_t fDistanceFromGoodClusterAndTrackAtLastPlane;
   Double_t fDistanceFromBestClusterAndTrackAtLastPlane;
 
-  TClonesArray *fMFTClusterArray[fMaxNPlanesMFT];         // array of clusters for the planes of the MFT
-  TClonesArray *fMFTClusterArrayFront[fMaxNPlanesMFT];    // array of front clusters for the planes of the MFT
-  TClonesArray *fMFTClusterArrayBack[fMaxNPlanesMFT];     // array of back clusters for the planes of the MFT
+  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
 
   Double_t fRAbsorberCut;  // in cm, corresponds to the radial position of a 3 degrees track at the end of the absorber (-503 cm)
   Double_t fLowPtCut;      // in GeV/c, the lower limit for the pt of a track in the muon spectrometer
   Int_t fNPlanesMFT;              // number of planes of the Vertex Telescope (Muon Internal Tracker) -> This should be taken from the new version of  AliVZERO2
   Int_t fNPlanesMFTAnalyzed;
   Int_t fNMaxMissingMFTClusters;  // max. number of MFT clusters which can be missed in the global fit procedure
-  Bool_t fIsPlaneMandatory[fMaxNPlanesMFT];    // specifies which MFT planes cannot be missed in the global fit procedure
+  Bool_t fIsPlaneMandatory[fNMaxPlanes];      // specifies which MFT planes cannot be missed in the global fit procedure
 
   Int_t fEv;               // current event being analyzed
   Int_t fLabelMC;          // MC label of the muon track reconstructed in the spectrometer
 
   Bool_t fIsClusterCompatible[10];       // here the clusters in the Muon Spectrometer are concerned
 
-  Double_t fZPlane[fMaxNPlanesMFT];        // z-position of the MFT planes (center of the support)
-  Double_t fRPlaneMax[fMaxNPlanesMFT];     // max radius of the MFT planes (the support)
-  Double_t fRPlaneMin[fMaxNPlanesMFT];     // min radius of the MFT planes (the support)
-
-  TH1D *fHistPtSpectrometer, *fHistPtMuonTrackWithGoodMatch, *fHistPtMuonTrackWithBadMatch;
-  TH1D *fHistRadiusEndOfAbsorber, *fHistNTracksAfterExtrapolation[fMaxNPlanesMFT];
-  TH1D *fHistNGoodClustersForFinalTracks, *fHistResearchRadius[fMaxNPlanesMFT];
-  TH1D *fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane;
-  TH1D *fHistDistanceGoodClusterFromTrackAtLastPlane;
-  TH1D *fHistChi2Cluster_GoodCluster[fMaxNPlanesMFT], *fHistChi2Cluster_BadCluster[fMaxNPlanesMFT];
-  TH1D *fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[fMaxNPlanesMFT], *fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[fMaxNPlanesMFT];
-
-  TNtuple *fNtuFinalCandidates1, *fNtuFinalBestCandidates1;
-  TNtuple *fNtuFinalCandidates2, *fNtuFinalBestCandidates2;
-
-  TGraph *fGrMFTPlane[fMaxNPlanesMFT][4];
-  TEllipse *fCircleExt[fMaxNPlanesMFT], *fCircleInt[fMaxNPlanesMFT];
-  TCanvas *fCanvas;
-
-  TLatex  *fTxtMuonHistory, *fTxtTrackGoodClusters, *fTxtTrackChi2[fMaxNPlanesMFT], *fTxtTrackFinalChi2, *fTxtFinalCandidates, *fTxtDummy;
-  TLatex  *fTxtAllClust, *fTxtClustGoodChi2, *fTxtClustMC, *fTxtClustOfTrack; 
-  TMarker *fMrkAllClust, *fMrkClustGoodChi2, *fMrkClustMC, *fMrkClustOfTrack;
-
-  Int_t fCountRealTracksAnalyzed; 
+  Double_t fZPlane[fNMaxPlanes];        // z-position of the MFT planes (center of the support)
+  Double_t fRPlaneMax[fNMaxPlanes];     // max radius of the MFT planes (the support)
+  Double_t fRPlaneMin[fNMaxPlanes];     // min radius of the MFT planes (the support)
+
+  TH1D *fHistPtSpectrometer, *fHistPtMuonTrackWithGoodMatch, *fHistPtMuonTrackWithBadMatch;     //
+  TH1D *fHistRadiusEndOfAbsorber, *fHistNTracksAfterExtrapolation[fNMaxPlanes];                        //
+  TH1D *fHistNGoodClustersForFinalTracks, *fHistResearchRadius[fNMaxPlanes];                   //
+  TH1D *fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane;         //
+  TH1D *fHistDistanceGoodClusterFromTrackAtLastPlane;                                          //
+  TH1D *fHistChi2Cluster_GoodCluster[fNMaxPlanes], *fHistChi2Cluster_BadCluster[fNMaxPlanes];   //
+  TH1D *fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[fNMaxPlanes];                        //
+  TH1D *fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[fNMaxPlanes];                  //
+
+  TNtuple *fNtuFinalCandidates;
+  TNtuple *fNtuFinalBestCandidates;
+
+  TGraph *fGrMFTPlane[4][20];             //!
+  TEllipse *fCircleExt[fNMaxPlanes], *fCircleInt[fNMaxPlanes];  //!
+  TCanvas *fCanvas;                       //!
+
+  TLatex  *fTxtMuonHistory, *fTxtTrackGoodClusters, *fTxtTrackChi2[fNMaxPlanes];      //!
+  TLatex  *fTxtTrackFinalChi2, *fTxtTrackMomentum, *fTxtFinalCandidates, *fTxtDummy;  //!
+  TLatex  *fTxtAllClust, *fTxtClustGoodChi2, *fTxtClustMC, *fTxtClustOfTrack;        //!
+  TMarker *fMrkAllClust, *fMrkClustGoodChi2, *fMrkClustMC, *fMrkClustOfTrack;         //!
+  Int_t fCountRealTracksAnalyzed;
+  Int_t fMaxNTracksToBeAnalyzed;
   Int_t fCountRealTracksWithRefMC; 
   Int_t fCountRealTracksWithRefMC_andTrigger;
   Int_t fCountRealTracksWithRefMC_andTrigger_andGoodPt;
@@ -212,45 +218,46 @@ protected:
   Int_t fCountRealTracksAnalyzedOfEvent;
   Int_t fCountRealTracksAnalyzedWithFinalCandidates;
 
-  Int_t fNClustersGlobalTrack[fMaxNPlanesMFT], fNDFGlobalTrack[fMaxNPlanesMFT];
+  Int_t fNClustersGlobalTrack[fNMaxPlanes], fNDFGlobalTrack[fNMaxPlanes];
   
-  TFile *fFileCluster;
-  TFile *fFileESD;
-  TFile *fFile_gAlice;
+  TFile *fFileCluster;   //!
+  TFile *fFileESD;       //!
+  TFile *fFile_gAlice;   //!
 
-  AliRunLoader *fRunLoader;
-  AliLoader *fMFTLoader;
-  AliMUONRecoCheck *fMuonRecoCheck;
+  AliRunLoader *fRunLoader;           //!
+  AliLoader *fMFTLoader;              //!
+  AliMUONRecoCheck *fMuonRecoCheck;   //!
 
-  TTree *fMFTClusterTree;
+  TTree *fMFTClusterTree;  //!
 
-  AliMUONTrack *fMuonTrackReco;           // muon track being analyzed
-  AliMuonForwardTrack *fCurrentTrack;     // muon extrapolated track being tested
+  AliMUONTrack *fMuonTrackReco;                 //! muon track being analyzed
+  AliMuonForwardTrack *fCurrentTrack;           //! muon extrapolated track being tested
+  AliMuonForwardTrack *fFinalBestCandidate;     //! best final candidate (if any)
   Bool_t fIsCurrentMuonTrackable;
-  Bool_t fIsGoodClusterInPlane[fMaxNPlanesMFT];
+  Bool_t fIsGoodClusterInPlane[fNMaxPlanes];
 
-  TClonesArray *fCandidateTracks;   // array of track we are going to build (starting from fMuonTrackReco)
+  TClonesArray *fCandidateTracks;   //! array of track we are going to build (starting from fMuonTrackReco)
 
-  AliMUONVTrackStore *fTrackStore;      // list of reconstructed MUON tracks 
-  AliMUONVTrackStore *fTrackRefStore;   // list of reconstructible MUON tracks
+  AliMUONVTrackStore *fTrackStore;      //! list of reconstructed MUON tracks 
+  AliMUONVTrackStore *fTrackRefStore;   //! list of reconstructible MUON tracks
   
   TIterator *fNextTrack;   //! Iterator for reading the MUON tracks
   
-  AliStack *fStack;
+  AliStack *fStack;  //!
 
-  AliMFT *fMFT;
-  AliMFTSegmentation *fSegmentation;
+  AliMFT *fMFT;                        //!
+  AliMFTSegmentation *fSegmentation;   //!
 
-  TFile *fOutputTreeFile;
-  TTree *fOutputEventTree;
+  TFile *fOutputTreeFile, *fOutputQAFile;   //
+  TTree *fOutputEventTree;                  //!
 
-  TClonesArray *fMuonForwardTracks;       // array of AliMuonForwardTrack
+  TClonesArray *fMuonForwardTracks;       //! array of AliMuonForwardTrack
 
   Int_t fMatchingMode;
   Double_t fMinResearchRadiusAtLastPlane;
 
-  AliGRPObject *fGRPData;              // Data from the GRP/GRP/Data CDB folder
-  AliRunInfo *fRunInfo;
+  AliGRPObject *fGRPData;              //! Data from the GRP/GRP/Data CDB folder
+  AliRunInfo *fRunInfo;                //!
 
   ClassDef(AliMuonForwardTrackFinder, 1); 
 
index ac23f1c..474c52c 100644 (file)
@@ -38,7 +38,8 @@ ClassImp(AliMuonForwardTrackPair)
 
 AliMuonForwardTrackPair::AliMuonForwardTrackPair():
   TObject(),
-  fMuonForwardTracks(0)
+  fMuonForwardTracks(0),
+  fKinemMC(0,0,0,0)
 {
 
   // default constructor
@@ -51,7 +52,8 @@ AliMuonForwardTrackPair::AliMuonForwardTrackPair():
 
 AliMuonForwardTrackPair::AliMuonForwardTrackPair(AliMuonForwardTrack *track0, AliMuonForwardTrack *track1):
   TObject(),
-  fMuonForwardTracks(0)
+  fMuonForwardTracks(0),
+  fKinemMC(0,0,0,0)
 {
 
   fMuonForwardTracks = new TClonesArray("AliMuonForwardTrack", 2);
@@ -59,13 +61,16 @@ AliMuonForwardTrackPair::AliMuonForwardTrackPair(AliMuonForwardTrack *track0, Al
   new ((*fMuonForwardTracks)[0]) AliMuonForwardTrack(*track0);
   new ((*fMuonForwardTracks)[1]) AliMuonForwardTrack(*track1);
 
+  SetKinemMC();
+
 }
 
 //====================================================================================================================================================
 
 AliMuonForwardTrackPair::AliMuonForwardTrackPair(const AliMuonForwardTrackPair& trackPair): 
   TObject(trackPair),
-  fMuonForwardTracks(trackPair.fMuonForwardTracks)
+  fMuonForwardTracks(trackPair.fMuonForwardTracks),
+  fKinemMC(trackPair.fKinemMC)
 {
 
   // copy constructor
@@ -135,11 +140,21 @@ Double_t AliMuonForwardTrackPair::GetMass(Double_t z, Int_t nClusters) {
   AliMUONTrackParam *param0 = ((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetTrackParamAtMFTCluster(idCluster[0]);
   AliMUONTrackParam *param1 = ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetTrackParamAtMFTCluster(idCluster[1]);
 
+  AliDebug(2, Form("MFT before extrap: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)", 
+                  param0->Px(), param0->Py(), param0->Pz(), 
+                  param1->Px(), param1->Py(), param1->Pz()));
+
   if (TMath::Abs(z)<1e6) {
+    AliDebug(2, Form("Extrapolating 1st muon from z = %f to z = %f", param0->GetZ(), z));
     AliMUONTrackExtrap::ExtrapToZCov(param0, z);
+    AliDebug(2, Form("Extrapolating 2nd muon from z = %f to z = %f", param1->GetZ(), z));
     AliMUONTrackExtrap::ExtrapToZCov(param1, z);
   }
 
+  AliDebug(2, Form("MFT after extrap: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)", 
+                  param0->Px(), param0->Py(), param0->Pz(), 
+                  param1->Px(), param1->Py(), param1->Pz()));
+
   momentum[0] = (param0->P());
   momentum[1] = (param1->P());
 
@@ -149,14 +164,9 @@ Double_t AliMuonForwardTrackPair::GetMass(Double_t z, Int_t nClusters) {
 
   dimu.SetE(TMath::Sqrt(mMu*mMu + momentum[0]*momentum[0]) + TMath::Sqrt(mMu*mMu + momentum[1]*momentum[1]));
 
-  dimu.SetPx(((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetTrackParamAtMFTCluster(idCluster[0])->Px()+
-            ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetTrackParamAtMFTCluster(idCluster[1])->Px());
-
-  dimu.SetPy(((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetTrackParamAtMFTCluster(idCluster[0])->Py()+
-            ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetTrackParamAtMFTCluster(idCluster[1])->Py());
-
-  dimu.SetPz(((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetTrackParamAtMFTCluster(idCluster[0])->Pz()+
-            ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetTrackParamAtMFTCluster(idCluster[1])->Pz());
+  dimu.SetPx(param0->Px() + param1->Px());
+  dimu.SetPy(param0->Py() + param1->Py());
+  dimu.SetPz(param0->Pz() + param1->Pz());
 
   return dimu.M();
 
@@ -177,9 +187,19 @@ Double_t AliMuonForwardTrackPair::GetMassWithoutMFT(Double_t x, Double_t y, Doub
   AliMUONTrackParam *param0 = ((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetTrackParamAtMUONCluster(idCluster[0]);
   AliMUONTrackParam *param1 = ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetTrackParamAtMUONCluster(idCluster[1]);
 
+  AliDebug(2, Form("MUON before extrap: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)", 
+                  param0->Px(), param0->Py(), param0->Pz(), 
+                  param1->Px(), param1->Py(), param1->Pz()));
+
+  AliDebug(2, Form("Extrapolating 1st muon from z = %f to z = %f", param0->GetZ(), z));
   AliMUONTrackExtrap::ExtrapToVertex(param0, x, y, z, 0., 0.);   // this should reproduce what is done in AliMUONESDInterface::MUONToESD(...) 
+  AliDebug(2, Form("Extrapolating 2nd muon from z = %f to z = %f", param1->GetZ(), z));
   AliMUONTrackExtrap::ExtrapToVertex(param1, x, y, z, 0., 0.);   // this should reproduce what is done in AliMUONESDInterface::MUONToESD(...) 
 
+  AliDebug(2, Form("MUON after extrap: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)", 
+                  param0->Px(), param0->Py(), param0->Pz(), 
+                  param1->Px(), param1->Py(), param1->Pz()));
+
   Double_t momentum[2] = {0}; 
 
   momentum[0] = param0->P();
@@ -201,32 +221,58 @@ Double_t AliMuonForwardTrackPair::GetMassWithoutMFT(Double_t x, Double_t y, Doub
 
 //====================================================================================================================================================
 
-Double_t AliMuonForwardTrackPair::GetMassMC() {
-
-  TLorentzVector dimu;
-
-  dimu.SetE(((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetMCTrackRef()->Energy() +
-           ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetMCTrackRef()->Energy());
-
-  dimu.SetPx(((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetMCTrackRef()->Px() +
-            ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetMCTrackRef()->Px());
-
-  dimu.SetPy(((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetMCTrackRef()->Py() +
-            ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetMCTrackRef()->Py());
+void AliMuonForwardTrackPair::SetKinemMC() {
 
-  dimu.SetPz(((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetMCTrackRef()->Pz() +
-            ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetMCTrackRef()->Pz());
+  AliDebug(2, Form("MC: 1st muon = (%f, %f, %f) 2nd muon = (%f, %f, %f)", 
+                  ((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetMCTrackRef()->Px(),
+                  ((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetMCTrackRef()->Py(),
+                  ((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetMCTrackRef()->Pz(),
+                  ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetMCTrackRef()->Px(),
+                  ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetMCTrackRef()->Py(),
+                  ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetMCTrackRef()->Pz()));
 
-  return dimu.M();
+  fKinemMC.SetE(((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetMCTrackRef()->Energy() +
+               ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetMCTrackRef()->Energy());
+  
+  fKinemMC.SetPx(((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetMCTrackRef()->Px() +
+                ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetMCTrackRef()->Px());
+  
+  fKinemMC.SetPy(((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetMCTrackRef()->Py() +
+                ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetMCTrackRef()->Py());
+  
+  fKinemMC.SetPz(((AliMuonForwardTrack*) fMuonForwardTracks->At(0))->GetMCTrackRef()->Pz() +
+                ((AliMuonForwardTrack*) fMuonForwardTracks->At(1))->GetMCTrackRef()->Pz());
 
 }
 
 //====================================================================================================================================================
 
+Bool_t AliMuonForwardTrackPair::IsResonance() {
+
+  Bool_t result = kFALSE;
 
+  Int_t labelMC[2] = {0};
+  Int_t codePDG[2] = {0};
+  
+  for (Int_t iTrack=0; iTrack<2; iTrack++) {
+    labelMC[iTrack] = ((AliMuonForwardTrack*) fMuonForwardTracks->At(iTrack))->GetParentMCLabel(0);
+    codePDG[iTrack] = ((AliMuonForwardTrack*) fMuonForwardTracks->At(iTrack))->GetParentPDGCode(0);
+  }
 
+  AliDebug(1, Form("Muons' mothers: (%d, %d)", labelMC[0], labelMC[1]));
 
+  if (labelMC[0]==labelMC[1] && codePDG[0]==codePDG[1] && (codePDG[0]==   113 ||
+                                                          codePDG[0]==   223 ||
+                                                          codePDG[0]==   333 ||
+                                                          codePDG[0]==   443 ||
+                                                          codePDG[0]==100443 ||
+                                                          codePDG[0]==   553 ||
+                                                          codePDG[0]==100553 ) ) result = kTRUE;
 
+  if (result) AliDebug(1, Form("Pair is a resonance %d", codePDG[0]));
 
+  return result; 
 
+}
 
+//====================================================================================================================================================
index 34e85b6..f922c3e 100644 (file)
@@ -36,16 +36,26 @@ public:
   virtual ~AliMuonForwardTrackPair() {}
 
   void SetTrack(Int_t iTrack, AliMuonForwardTrack *track);
-  AliMuonForwardTrack* GetTrack(Int_t iTrack) { if (iTrack==0 || iTrack==1) return (AliMuonForwardTrack*) fMuonForwardTracks->At(iTrack); else return NULL; }
+  AliMuonForwardTrack* GetTrack(Int_t iTrack) { 
+    if (iTrack==0 || iTrack==1) return (AliMuonForwardTrack*) fMuonForwardTracks->At(iTrack); 
+    else return NULL; 
+  }
+
+  void SetKinemMC();
 
   Double_t GetWeightedOffset(Double_t x, Double_t y, Double_t z);
   Double_t GetMass(Double_t z, Int_t nClusters=-1);
   Double_t GetMassWithoutMFT(Double_t x, Double_t y, Double_t z, Int_t nClusters=-1);
-  Double_t GetMassMC();
+  Double_t GetMassMC()     { return fKinemMC.M(); }
+  Double_t GetRapidityMC() { return fKinemMC.Rapidity(); }
+  Double_t GetPtMC()       { return fKinemMC.Pt(); }
+
+  Bool_t IsResonance();
 
 protected:
 
   TClonesArray *fMuonForwardTracks;
+  TLorentzVector fKinemMC;
 
   ClassDef(AliMuonForwardTrackPair,1)
     
index 9eec69d..b429f3c 100644 (file)
@@ -25,7 +25,7 @@
 # SHLIBS - Shared Libraries and objects for linking (Executables only)           #
 #--------------------------------------------------------------------------------#
 
-set ( SRCS  AliMFTPlane.cxx AliMFTSegmentation.cxx AliMFTDigit.cxx AliMFTCluster.cxx )
+set ( SRCS  AliMFTConstants.cxx AliMFTPlane.cxx AliMFTSegmentation.cxx AliMFTDigit.cxx AliMFTCluster.cxx )
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
 
index 0ebfe5c..3f4362f 100644 (file)
@@ -6,6 +6,7 @@
 #pragma link off all classes;
 #pragma link off all functions;
  
+#pragma link C++ class  AliMFTConstants+;
 #pragma link C++ class  AliMFTPlane+;
 #pragma link C++ class  AliMFTSegmentation+;
 #pragma link C++ class  AliMFTDigit+;
index 4a3528a..a86e6f9 100644 (file)
@@ -10,10 +10,10 @@ void SetMFTGeometry() {
 
   const Int_t nPlanes = 5;
   
-  const Float_t zCenter[nPlanes]          = {  -50.0,   -58.0,   -66.0,   -74.0,   -82.0 };   // expressed in cm
+  const Float_t zCenter[nPlanes]          = {  -50.0,   -58.0,   -66.0,   -72.0,   -76.0 };   // expressed in cm
                                          
-  const Float_t rMin[nPlanes]             = {   2.00,    2.31,    2.66,    3.01,    3.36 };   // expressed in cm  
-  const Float_t rMax[nPlanes]             = {   9.70,   11.11,   12.52,   13.93,   15.34 };   // expressed in cm
+  const Float_t rMin[nPlanes]             = {   2.20,    2.60,    3.00,    3.30,    3.60 };   // expressed in cm  
+  const Float_t rMax[nPlanes]             = {   9.70,   11.00,   12.40,   13.40,   14.00 };   // expressed in cm
                                          
   const Float_t pixelSizeX[nPlanes]       = { 20.e-4,  20.e-4,  20.e-4,  20.e-4,  20.e-4 };   // expressed in cm
   const Float_t pixelSizeY[nPlanes]       = { 20.e-4,  20.e-4,  20.e-4,  20.e-4,  20.e-4 };   // expressed in cm