Changes to EventReconstructor...:
authorgosset <gosset@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 30 Jun 2000 10:15:48 +0000 (10:15 +0000)
committergosset <gosset@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 30 Jun 2000 10:15:48 +0000 (10:15 +0000)
precision fit with multiple Coulomb scattering;
extrapolation to vertex with Branson correction in absorber (JPC)

MUON/AliMUONEventReconstructor.cxx
MUON/AliMUONEventReconstructor.h
MUON/AliMUONSegment.cxx
MUON/AliMUONSegment.h
MUON/AliMUONTrack.cxx
MUON/AliMUONTrack.h
MUON/AliMUONTrackParam.cxx
MUON/AliMUONTrackParam.h

index 3c0557ff6733800d02ccac8c4183d7141d8d6c66..714200fe28eba6e464bc113c50694667ea6e0a06 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.4  2000/06/27 14:11:36  gosset
+Corrections against violations of coding conventions
+
 Revision 1.3  2000/06/16 07:27:08  gosset
 To remove problem in running RuleChecker, like in MUON-dev
 
@@ -67,6 +70,7 @@ Addition of files for track reconstruction in C++
 #include "AliMUONChamber.h"
 #include "AliMUONTrackHit.h"
 #include "AliRun.h"
+#include "TParticle.h"
 
 #ifndef WIN32 
 # define initfield initfield_
@@ -849,6 +853,7 @@ Int_t AliMUONEventReconstructor::MakeTrackCandidatesWithTwoSegments(AliMUONSegme
     // pointer to segment
     endSegment = (AliMUONSegment*) ((*fSegmentsPtr[endStation])[iEndSegment]);
     // test compatibility between current segment and "extrapSegment"
+    // 4 because 4 quantities in chi2
     if ((endSegment->
         NormalizedChi2WithSegment(extrapSegment,
                                   fMaxSigma2Distance)) <= 4.0) {
@@ -909,6 +914,7 @@ Int_t AliMUONEventReconstructor::MakeTrackCandidatesWithOneSegmentAndOnePoint(Al
       // pointer to HitForRec
       hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
       // test compatibility between current HitForRec and "extrapHitForRec"
+      // 2 because 2 quantities in chi2
       if ((hit->
           NormalizedChi2WithHitForRec(extrapHitForRec,
                                       fMaxSigma2Distance)) <= 2.0) {
@@ -973,23 +979,29 @@ void AliMUONEventReconstructor::MakeTrackCandidates(void)
 void AliMUONEventReconstructor::FollowTracks(void)
 {
   // Follow tracks in stations(1..) 3, 2 and 1
+  // too long: should be made more modular !!!!
   AliMUONHitForRec *bestHit, *extrapHit, *hit;
   AliMUONSegment *bestSegment, *extrapSegment, *segment;
-  AliMUONTrack *track;
-  AliMUONTrackParam *trackParam1, trackParam[2];
+  AliMUONTrack *track, *nextTrack;
+  AliMUONTrackParam *trackParam1, trackParam[2], trackParamVertex;
   Int_t ch, chBestHit, iHit, iSegment, station, trackIndex;
-  Double_t bestChi2, chi2, dZ1, dZ2, maxSigma2Distance, mcsFactor;
+  Double_t bestChi2, chi2, dZ1, dZ2, dZ3, maxSigma2Distance, mcsFactor;
   AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
-  maxSigma2Distance = 4. * fMaxSigma2Distance; // sigma2cut increased by 4 !!!!
+  // local maxSigma2Distance, for easy increase in testing
+  maxSigma2Distance = fMaxSigma2Distance;
   if (fPrintLevel >= 2)
     cout << "enter FollowTracks" << endl;
   // Loop over track candidates
-  for (trackIndex = 0; trackIndex < fNRecTracks; trackIndex++) {
+  track = (AliMUONTrack*) fRecTracksPtr->First();
+  trackIndex = -1;
+  while (track) {
+    // Follow function for each track candidate ????
+    trackIndex++;
+    nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track
     if (fPrintLevel >= 2)
       cout << "FollowTracks: track candidate(0..): " << trackIndex << endl;
-    // function for each track candidate ????
-    track = (AliMUONTrack*) ((*fRecTracksPtr)[trackIndex]);
     // Fit track candidate from vertex at X = Y = 0 
+    track->SetFitMCS(0);  // fit without Multiple Scattering
     track->Fit(track->GetTrackParamAtVertex(), 3);
     if (fPrintLevel >= 10) {
       cout << "FollowTracks: track candidate(0..): " << trackIndex
@@ -997,7 +1009,8 @@ void AliMUONEventReconstructor::FollowTracks(void)
       track->RecursiveDump();
     }
     // Loop over stations(1..) 3, 2 and 1
-    // something SPECIAL for stations 2 and 1 for majority coincidence ????
+    // something SPECIAL for stations 2 and 1 for majority 3 coincidence ????
+    // otherwise: majority coincidence 2 !!!!
     for (station = 2; station >= 0; station--) {
       // Track parameters at first track hit (smallest Z)
       trackParam1 = ((AliMUONTrackHit*)
@@ -1015,9 +1028,15 @@ void AliMUONEventReconstructor::FollowTracks(void)
       // Z difference between the two previous stations
       dZ2 = (&(pMUON->Chamber(2 * station + 2)))->Z() -
        (&(pMUON->Chamber(2 * station + 4)))->Z();
-      extrapSegment->SetBendingCoorReso2(fBendingResolution);
-      extrapSegment->SetNonBendingCoorReso2(fNonBendingResolution);
-      extrapSegment->UpdateFromStationTrackParam(trackParam, mcsFactor, dZ1, dZ2);
+      // Z difference between the two chambers in the previous station
+      dZ3 = (&(pMUON->Chamber(2 * station)))->Z() -
+       (&(pMUON->Chamber(2 * station + 1)))->Z();
+      extrapSegment->SetBendingCoorReso2(fBendingResolution * fBendingResolution);
+      extrapSegment->
+       SetNonBendingCoorReso2(fNonBendingResolution * fNonBendingResolution);
+      extrapSegment->UpdateFromStationTrackParam
+       (trackParam, mcsFactor, dZ1, dZ2, dZ3, station,
+        trackParam1->GetInverseBendingMomentum());
       bestChi2 = 5.0;
       bestSegment = NULL;
       if (fPrintLevel >= 10) {
@@ -1100,26 +1119,125 @@ void AliMUONEventReconstructor::FollowTracks(void)
          }
        }
        else {
-         fRecTracksPtr->RemoveAt(trackIndex); // Remove candidate
+         // Remove current track candidate and update fNRecTracks
+         // To be checked: recursive delete of TrackHit's !!!!
+         // For cleaner implementation: call track->Remove()
+         // to be coded, with all cleanings,
+         // including links between HitForRec's and TrackHit's !!!!
+         fRecTracksPtr->Remove(track);
+         fNRecTracks--;
+         delete extrapSegment;
          break; // stop the search for this candidate:
          // exit from the loop over station
        }
       }
+      delete extrapSegment;
       // Sort track hits according to increasing Z
       track->GetTrackHitsPtr()->Sort();
       // Update track parameters at first track hit (smallest Z)
       trackParam1 = ((AliMUONTrackHit*)
                     (track->GetTrackHitsPtr()->First()))->GetTrackParam();
-      // Track fit from first track hit varying X and Y
+      // Track fit from first track hit varying X and Y,
+      // with multiple Coulomb scattering if all stations
+      if (station == 0) track->SetFitMCS(1);
+      else track->SetFitMCS(0);
       track->Fit(trackParam1, 5);
       if (fPrintLevel >= 10) {
        cout << "FollowTracks: track candidate(0..): " << trackIndex
             << " after fit from station(0..): " << station << " to 4" << endl;
        track->RecursiveDump();
       }
-      delete extrapSegment;
+      // Track extrapolation to the vertex through the absorber (Branson)
+      if (station == 0) {
+       trackParamVertex = *trackParam1;
+       (&trackParamVertex)->ExtrapToVertex();
+       track->SetTrackParamAtVertex(&trackParamVertex);
+      }
+      if (fPrintLevel >= 1) {
+       cout << "FollowTracks: track candidate(0..): " << trackIndex
+            << " after fit from station(0..): " << station << " to 4" << endl;
+       //      track->RecursiveDump(); // CCC
+      }
     } // for (station = 2;...
-  } // for (trackIndex = 0;...
+    // go really to next track
+    track = nextTrack;
+  } // while (track)
+  // Compression of track array (necessary after Remove ????)
+  fRecTracksPtr->Compress();
+  return;
+}
+
+  //__________________________________________________________________________
+void AliMUONEventReconstructor::EventDump(void)
+{
+  // Dump reconstructed event (track parameters at vertex and at first hit),
+  // and the particle parameters
+
+  AliMUONTrack *track;
+  AliMUONTrackParam *trackParam, *trackParam1;
+  TClonesArray *particles; // pointer to the particle list
+  TParticle *p;
+  Double_t bendingSlope, nonBendingSlope, pYZ;
+  Double_t pX, pY, pZ, x, y, z, c;
+  Int_t np, trackIndex, nTrackHits;
+  if (fPrintLevel >= 1) cout << "****** enter EventDump ******" << endl;
+  if (fPrintLevel >= 1) {
+    cout << " Number of Reconstructed tracks :" <<  fNRecTracks << endl; 
+  }
+  fRecTracksPtr->Compress(); // for simple loop without "Next" since no hole
+  // Loop over reconstructed tracks
+  for (trackIndex = 0; trackIndex < fNRecTracks; trackIndex++) {
+    if (fPrintLevel >= 1)
+      cout << " track number: " << trackIndex << endl;
+    // function for each track for modularity ????
+    track = (AliMUONTrack*) ((*fRecTracksPtr)[trackIndex]);
+    nTrackHits = track->GetNTrackHits();
+    if (fPrintLevel >= 1)
+      cout << " number of track hits: " << nTrackHits << endl;
+    // track parameters at Vertex
+    trackParam = track->GetTrackParamAtVertex();
+    x = trackParam->GetNonBendingCoor();
+    y = trackParam->GetBendingCoor();
+    z = trackParam->GetZ();
+    bendingSlope = trackParam->GetBendingSlope();
+    nonBendingSlope = trackParam->GetNonBendingSlope();
+    pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
+    pZ = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
+    pX = pZ * nonBendingSlope;
+    pY = pZ * bendingSlope;
+    c = TMath::Sign(1.0, trackParam->GetInverseBendingMomentum());
+    if (fPrintLevel >= 1)
+      printf(" track parameters at Vertex z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
+            z, x, y, pX, pY, pZ, c);
+
+    // track parameters at first hit
+    trackParam1 = ((AliMUONTrackHit*)
+                  (track->GetTrackHitsPtr()->First()))->GetTrackParam();
+    x = trackParam1->GetNonBendingCoor();
+    y = trackParam1->GetBendingCoor();
+    z = trackParam1->GetZ();
+    bendingSlope = trackParam1->GetBendingSlope();
+    nonBendingSlope = trackParam1->GetNonBendingSlope();
+    pYZ = 1/TMath::Abs(trackParam1->GetInverseBendingMomentum());
+    pZ = pYZ/TMath::Sqrt(1.0 + bendingSlope * bendingSlope);
+    pX = pZ * nonBendingSlope;
+    pY = pZ * bendingSlope;
+    c = TMath::Sign(1.0, trackParam1->GetInverseBendingMomentum());
+    if (fPrintLevel >= 1)
+      printf(" track parameters at z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
+            z, x, y, pX, pY, pZ, c);
+  }
+  // informations about generated particles
+  particles = gAlice->Particles();
+  np = particles->GetEntriesFast();
+  printf(" **** number of generated particles: %d  \n", np);
+  
+  for (Int_t iPart = 0; iPart < np; iPart++) {
+    p = (TParticle*) particles->UncheckedAt(iPart);
+    printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n",
+          iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode());    
+  }
   return;
 }
 
index 396821f4df358ef0251a89fc0dc58d440b84fe0f..bddc9521a0707cf91c38f9338724b17c2bf48132 100644 (file)
@@ -64,6 +64,7 @@ class AliMUONEventReconstructor : public TObject {
   Double_t GetImpactParamFromBendingMomentum(Double_t BendingMomentum);
   Double_t GetBendingMomentumFromImpactParam(Double_t ImpactParam);
   void EventReconstruct(void);
+  void EventDump(void);  // dump reconstructed event
 
  protected:
 
index 36406f7b6aab9a7b56d55e2fa0bbd54046cda5e0..01671570e80cc2e02eb63c71086850806c7e275b 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.3  2000/06/25 13:06:39  hristov
+Inline functions moved from *.cxx to *.h files instead of forward declarations
+
 Revision 1.2  2000/06/15 07:58:48  morsch
 Code from MUON-dev joined
 
@@ -229,7 +232,7 @@ AliMUONHitForRec* AliMUONSegment::CreateHitForRecFromLinearExtrapToChamber (Int_
 }
 
   //__________________________________________________________________________
-void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2)
+void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2, Double_t Dz3, Int_t Station, Double_t InverseMomentum)
 {
   // Fill data members with values calculated from the array of track parameters
   // pointed to by "TrackParam" (index = 0 and 1 for first and second chambers
@@ -239,28 +242,108 @@ void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam,
   // with one chamber for the coordinate, two chambers for the angle,
   // due to the arrangement in stations.
   // Resolution coming from:
-  // coordinate in closest station at "Dz1",
+  // coordinate in closest station at "Dz1" from current "Station",
   // slope between closest stations, with "Dz2" interval between them,
-  // extrapolation over "Dz" from closest station.
+  // interval "Dz3" between chambers of closest station,
+  // extrapolation over "Dz1" from closest station,
+  // "InverseMomentum".
   // When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
   // are assumed to be filled
   // with the variance on bending and non bending coordinates.
+  // The "road" is parametrized from the old reco_muon.F
+  // with 8 cm between stations.
   AliMUONTrackParam *param0;
   Double_t cReso2, sReso2;
+  // parameters to define the widths of the searching roads in station 0,1,2
+  // width = p0 + p1/ (momentum)^2
+  //                  station number:        0         1          2
+  static Double_t p0BendingCoor[3] =     { 6.43e-2, 1.64e-2,   0.034 };   
+  static Double_t p1BendingCoor[3] =     {    986.,    821.,    446. };  
+  static Double_t p0BendingSlope[3] =    { 3.54e-6, 3.63e-6,  3.6e-6 };  
+  static Double_t p1BendingSlope[3] =    { 4.49e-3,  4.8e-3,   0.011 };  
+  static Double_t p0NonBendingCoor[3] =  { 4.66e-2, 4.83e-2,   0.049 };   
+  static Double_t p1NonBendingCoor[3] =  {   1444.,    866.,    354. };  
+  static Double_t p0NonBendingSlope[3] = { 6.14e-4, 6.49e-4, 6.85e-4 };  
+  static Double_t p1NonBendingSlope[3] = {      0.,      0.,      0. };  
   param0 = &(TrackParam[0]);
+
+// OLD version
+//   // Bending plane
+//   fBendingCoor = param0->GetBendingCoor(); // coordinate
+//   fBendingSlope = param0->GetBendingSlope(); // slope
+//   cReso2 = fBendingCoorReso2;
+//   sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
+//   fBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
+//   fBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
+//   // Non bending plane
+//   fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
+//   fNonBendingSlope = param0->GetNonBendingSlope(); // slope
+//   cReso2 = fNonBendingCoorReso2;
+//   sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
+//   fNonBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
+//   fNonBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
+
+  // Coordinate and slope
   // Bending plane
   fBendingCoor = param0->GetBendingCoor(); // coordinate
   fBendingSlope = param0->GetBendingSlope(); // slope
-  cReso2 = fBendingCoorReso2;
-  sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
-  fBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
-  fBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
   // Non bending plane
   fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
   fNonBendingSlope = param0->GetNonBendingSlope(); // slope
+
+  // Resolutions
+  // cReso2 and sReso2 have to be subtracted here from the parametrization
+  // because they are added in the functions "NormalizedChi2WithSegment"
+  // and "NormalizedChi2WithHitForRec"
+  // Bending plane
+  cReso2 = fBendingCoorReso2;
+  sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
+  fBendingCoorReso2 = p0BendingCoor[Station] + p1BendingCoor[Station]*InverseMomentum*InverseMomentum - cReso2;
+  fBendingSlopeReso2 = p0BendingSlope[Station] + p1BendingSlope[Station]*InverseMomentum*InverseMomentum - sReso2;
+  // Non bending plane
   cReso2 = fNonBendingCoorReso2;
-  sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
-  fNonBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
-  fNonBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
+  sReso2 =  (2. * cReso2 )/ (Dz3*Dz3) ;
+  fNonBendingCoorReso2 = p0NonBendingCoor[Station] + p1NonBendingCoor[Station]*InverseMomentum*InverseMomentum - cReso2;
+  fNonBendingSlopeReso2 = p0NonBendingSlope[Station] + p1NonBendingSlope[Station]*InverseMomentum*InverseMomentum - sReso2;
   return;
 }
+
+// OLD function, with roads automatically calculated instead from being parametrized
+// kept because it would be a better solution,
+// if one can really find the right values.
+//   //__________________________________________________________________________
+// void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2)
+// {
+//   // Fill data members with values calculated from the array of track parameters
+//   // pointed to by "TrackParam" (index = 0 and 1 for first and second chambers
+//   // of the station, respectively).
+//   // Multiple Coulomb scattering is taking into account with "MCSfactor"
+//   // corresponding to one chamber,
+//   // with one chamber for the coordinate, two chambers for the angle,
+//   // due to the arrangement in stations.
+//   // Resolution coming from:
+//   // coordinate in closest station at "Dz1",
+//   // slope between closest stations, with "Dz2" interval between them,
+//   // extrapolation over "Dz" from closest station.
+//   // When called, "fBendingCoorReso2" and "fNonBendingCoorReso2"
+//   // are assumed to be filled
+//   // with the variance on bending and non bending coordinates.
+//   AliMUONTrackParam *param0;
+//   Double_t cReso2, sReso2;
+//   param0 = &(TrackParam[0]);
+//   // Bending plane
+//   fBendingCoor = param0->GetBendingCoor(); // coordinate
+//   fBendingSlope = param0->GetBendingSlope(); // slope
+//   cReso2 = fBendingCoorReso2;
+//   sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
+//   fBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
+//   fBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
+//   // Non bending plane
+//   fNonBendingCoor = param0->GetNonBendingCoor(); // coordinate
+//   fNonBendingSlope = param0->GetNonBendingSlope(); // slope
+//   cReso2 = fNonBendingCoorReso2;
+//   sReso2 = 2.0 * cReso2 / Dz2 / Dz2;
+//   fNonBendingCoorReso2 = cReso2 + (sReso2 + MCSfactor) * Dz1 * Dz1;
+//   fNonBendingSlopeReso2 = sReso2 + 2.0 * MCSfactor;
+//   return;
+// }
index 522f2034d0a85f54fd97faacaec7fb695c947c4b..5ce4e10b4b75408c9d912e4c7a2ad8f0888a666a 100644 (file)
@@ -54,7 +54,7 @@ class AliMUONSegment : public TObject {
   AliMUONSegment* CreateSegmentFromLinearExtrapToStation (Int_t Station, Double_t MCSfactor);
   Double_t NormalizedChi2WithSegment(AliMUONSegment* Segment, Double_t Sigma2Cut);
   AliMUONHitForRec* CreateHitForRecFromLinearExtrapToChamber (Int_t Chamber, Double_t MCSfactor);
-  void UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2);
+  void UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor, Double_t Dz1, Double_t Dz2, Double_t Dz3, Int_t Station, Double_t InverseMomentum);
 
   // What is necessary for sorting TClonesArray's; sufficient too ????
   Bool_t IsSortable() const { return kTRUE; }
index ba367af8bc63cd2c7b9d87a2d1e02220f2fc9aca..d1f77f1b18f68ddac65c6ca678c4350188d8ab26 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.3  2000/06/25 13:23:28  hristov
+stdlib.h needed for non-Linux compilation
+
 Revision 1.2  2000/06/15 07:58:48  morsch
 Code from MUON-dev joined
 
@@ -40,6 +43,8 @@ Addition of files for track reconstruction in C++
 
 #include <TClonesArray.h>
 #include <TMinuit.h>
+#include <TMath.h>
+#include <TMatrix.h>
 
 #include "AliMUONEventReconstructor.h" 
 #include "AliMUONHitForRec.h" 
@@ -48,11 +53,15 @@ Addition of files for track reconstruction in C++
 
 #include <stdlib.h>
 
+// variables to be known from minimization functions
 static AliMUONTrack *trackBeingFitted;
 static AliMUONTrackParam *trackParamBeingFitted;
 
-// Function to be minimized with Minuit
+// Functions to be minimized with Minuit
 void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
+void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
+
+Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit);
 
 ClassImp(AliMUONTrack) // Class implementation in ROOT context
 
@@ -68,6 +77,7 @@ AliMUONTrack::AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegmen
   AddSegment(EndSegment); // add hits from EndSegment
   fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
   SetTrackParamAtVertex(); // set track parameters at vertex
+  fFitMCS = 0;
   return;
 }
 
@@ -83,6 +93,7 @@ AliMUONTrack::AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec,
   AddHitForRec(HitForRec); // add HitForRec
   fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
   SetTrackParamAtVertex(); // set track parameters at vertex
+  fFitMCS = 0;
   return;
 }
 
@@ -97,6 +108,16 @@ AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& MUONTrack)
     return *this;
 }
 
+  //__________________________________________________________________________
+void AliMUONTrack::SetFitMCS(Int_t FitMCS)
+{
+  // Set track fit option with or without multiple Coulomb scattering
+  // from "FitMCS" argument: 0 without, 1 with
+  fFitMCS = FitMCS;
+  // better implementation with enum(with, without) ????
+  return;
+}
+
 // Inline functions for Get and Set: inline removed because it does not work !!!!
 AliMUONTrackParam* AliMUONTrack::GetTrackParamAtVertex(void) {
   // Get pointer to fTrackParamAtVertex
@@ -152,22 +173,26 @@ void AliMUONTrack::Fit(AliMUONTrackParam *TrackParam, Int_t NParam)
   TMinuit *minuit = new TMinuit(5);
   trackBeingFitted = this; // for the track to be known from the function to minimize
   trackParamBeingFitted = TrackParam; // for the track parameters to be known from the function to minimize; possible to use only Minuit parameters ????
+  // try to use TVirtualFitter to get this feature automatically !!!!
   minuit->mninit(5, 10, 7); // sysrd, syswr, syssa: useful ????
   // how to go faster ???? choice of Minuit parameters like EDM ????
-  minuit->SetFCN(TrackChi2);
+  // choice of function to be minimized according to fFitMCS
+  if (fFitMCS == 0) minuit->SetFCN(TrackChi2);
+  else minuit->SetFCN(TrackChi2MCS);
   minuit->SetPrintLevel(1); // More printing !!!!
-  // set first 3 parameters (try with no limits for the search: min=max=0)
+  // set first 3 parameters
+  // could be tried with no limits for the search (min=max=0) ????
   minuit->mnparm(0, "InvBenP",
                 TrackParam->GetInverseBendingMomentum(),
-                0.003, 0.0, 0.0, error);
+                0.003, -0.4, 0.4, error);
   minuit->mnparm(1, "BenS",
                 TrackParam->GetBendingSlope(),
-                0.001, 0.0, 0.0, error);
+                0.001, -0.5, 0.5, error);
   minuit->mnparm(2, "NonBenS",
                 TrackParam->GetNonBendingSlope(),
-                0.001, 0.0, 0.0, error);
+                0.001, -0.5, 0.5, error);
   if (NParam == 5) {
-    // set last 2 parameters (try with no limits for the search: min=max=0)
+    // set last 2 parameters (no limits for the search: min=max=0)
     minuit->mnparm(3, "X",
                   TrackParam->GetNonBendingCoor(),
                   0.03, 0.0, 0.0, error);
@@ -277,6 +302,7 @@ void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Para
   // and their current values in array pointed to by "Param".
   // Assumes that the track hits are sorted according to increasing Z.
   // Track parameters at each TrackHit are updated accordingly.
+  // Multiple Coulomb scattering is not taken into account
   AliMUONTrackHit* hit;
   AliMUONTrackParam param1;
   Int_t hitNumber;
@@ -314,3 +340,151 @@ void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Para
       dY * dY / hit->GetHitForRecPtr()->GetBendingReso2();
   }
 }
+
+  //__________________________________________________________________________
+void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag)
+{
+  // Return the "Chi2" to be minimized with Minuit for track fitting,
+  // with "NParam" parameters
+  // and their current values in array pointed to by "Param".
+  // Assumes that the track hits are sorted according to increasing Z.
+  // Track parameters at each TrackHit are updated accordingly.
+  // Multiple Coulomb scattering is taken into account with covariance matrix.
+  AliMUONTrackParam param1;
+  Chi2 = 0.0; // initialize Chi2
+  // copy of track parameters to be fitted
+  param1 = *trackParamBeingFitted;
+  // Minuit parameters to be fitted into this copy
+  param1.SetInverseBendingMomentum(Param[0]);
+  param1.SetBendingSlope(Param[1]);
+  param1.SetNonBendingSlope(Param[2]);
+  if (NParam == 5) {
+    param1.SetNonBendingCoor(Param[3]);
+    param1.SetBendingCoor(Param[4]);
+  }
+
+  AliMUONTrackHit* hit, hit1, hit2, hit3;
+  Bool_t GoodDeterminant;
+  Int_t hitNumber, hitNumber1, hitNumber2, hitNumber3;
+  Double_t zHit[10], paramBendingCoor[10], paramNonBendingCoor[10], ap[10];
+  Double_t hitBendingCoor[10], hitNonBendingCoor[10];
+  Double_t hitBendingReso2[10], hitNonBendingReso2[10];
+  Int_t numberOfHit = TMath::Min(trackBeingFitted->GetNTrackHits(), 10);
+  TMatrix *covBending = new TMatrix(numberOfHit, numberOfHit);
+  TMatrix *covNonBending = new TMatrix(numberOfHit, numberOfHit);
+
+  // Predicted coordinates and  multiple scattering angles are first calculated
+  for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) {
+    hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber];
+    zHit[hitNumber] = hit->GetHitForRecPtr()->GetZ();
+    // do something special if 2 hits with same Z ????
+    // security against infinite loop ????
+    (&param1)->ExtrapToZ(zHit[hitNumber]); // extrapolation
+    hit->SetTrackParam(&param1);
+    paramBendingCoor[hitNumber]= (&param1)->GetBendingCoor();
+    paramNonBendingCoor[hitNumber]= (&param1)->GetNonBendingCoor();
+    hitBendingCoor[hitNumber]= hit->GetHitForRecPtr()->GetBendingCoor();
+    hitNonBendingCoor[hitNumber]= hit->GetHitForRecPtr()->GetNonBendingCoor();
+    hitBendingReso2[hitNumber]= hit->GetHitForRecPtr()->GetBendingReso2();
+    hitNonBendingReso2[hitNumber]= hit->GetHitForRecPtr()->GetNonBendingReso2();
+    ap[hitNumber] = MultipleScatteringAngle2(hit); // multiple scatt. angle ^2  
+  }
+
+  // Calculates the covariance matrix
+  for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) {    
+    for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) {
+      (*covBending)(hitNumber1, hitNumber2) = 0;
+      (*covBending)(hitNumber2, hitNumber1) = 0;
+      if (hitNumber1 == hitNumber2){ // diagonal elements
+       (*covBending)(hitNumber2, hitNumber1) =
+         (*covBending)(hitNumber2, hitNumber1) + hitBendingReso2[hitNumber1];
+      }
+      // Multiple Scattering...  loop on upstream chambers ??
+      for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++){     
+       (*covBending)(hitNumber2, hitNumber1) =
+         (*covBending)(hitNumber2, hitNumber1) +
+         ((zHit[hitNumber1] - zHit[hitNumber3]) *
+          (zHit[hitNumber2] - zHit[hitNumber3]) * ap[hitNumber3]); 
+      }  
+      (*covNonBending)(hitNumber1, hitNumber2) = 0;
+      (*covNonBending)(hitNumber2, hitNumber1) =
+       (*covBending)(hitNumber2, hitNumber1);
+      if (hitNumber1 == hitNumber2) {  // diagonal elements
+       (*covNonBending)(hitNumber2, hitNumber1) =
+         (*covNonBending)(hitNumber2, hitNumber1) -
+         hitBendingReso2[hitNumber1] + hitNonBendingReso2[hitNumber1] ;
+      }      
+    }
+  }
+
+  // Inverts covariance matrix 
+  GoodDeterminant = kTRUE;
+  if (covBending->Determinant() != 0) {
+    covBending->Invert();
+  } else {
+    GoodDeterminant = kFALSE;
+    cout << "Warning in ChiMCS  Determinant Bending=0: " << endl;  
+  }
+  if (covNonBending->Determinant() != 0){
+    covNonBending->Invert();
+  } else {
+    GoodDeterminant = kFALSE;
+    cout << "Warning in ChiMCS  Determinant non Bending=0: " << endl;  
+  }
+  
+  // Calculates Chi2
+  if (GoodDeterminant) { // with Multiple Scattering if inversion correct
+    for (hitNumber1=0; hitNumber1 < numberOfHit ; hitNumber1++){ 
+      for (hitNumber2=0; hitNumber2 < numberOfHit; hitNumber2++){
+       Chi2 = Chi2 +
+         ((*covBending)(hitNumber2, hitNumber1) * 
+          (hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) *
+          (hitBendingCoor[hitNumber2] - paramBendingCoor[hitNumber2]));
+       Chi2 = Chi2 +
+         ((*covNonBending)(hitNumber2, hitNumber1) *
+          (hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) *
+          (hitNonBendingCoor[hitNumber2] - paramNonBendingCoor[hitNumber2]));
+      }
+    }
+  } else {  // without Multiple Scattering if inversion impossible
+    for (hitNumber1=0; hitNumber1 < numberOfHit ; hitNumber1++) { 
+      Chi2 = Chi2 +
+       ((hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) *
+        (hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) /
+        hitBendingReso2[hitNumber1]);
+      Chi2 = Chi2 +
+       ((hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) *
+        (hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) /
+        hitNonBendingReso2[hitNumber1]);      
+    }
+  }
+  
+  delete covBending;
+  delete covNonBending;
+}
+
+Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit)
+{
+  // Returns square of multiple Coulomb scattering angle
+  // at TrackHit pointed to by "TrackHit"
+  Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2;
+  Double_t varMultipleScatteringAngle;
+  AliMUONTrackParam *param = TrackHit->GetTrackParam();
+  slopeBending = param->GetBendingSlope();
+  slopeNonBending = param->GetNonBendingSlope();
+  // thickness in radiation length for the current track,
+  // taking local angle into account
+  radiationLength =
+    trackBeingFitted->GetEventReconstructor()->GetChamberThicknessInX0() *
+    TMath::Sqrt(1.0 +
+               slopeBending * slopeBending + slopeNonBending * slopeNonBending);
+  inverseBendingMomentum2 = 
+    param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum();
+  inverseTotalMomentum2 =
+    inverseBendingMomentum2 * (1.0 + slopeBending*slopeBending) /
+    (1.0 + slopeBending *slopeBending + slopeNonBending * slopeNonBending); 
+  varMultipleScatteringAngle = 0.0136 * (1.0 + 0.038 * TMath::Log(radiationLength));
+  varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength *
+    varMultipleScatteringAngle * varMultipleScatteringAngle;
+  return varMultipleScatteringAngle;
+}
index b528237fb18ec3afb316e6e2e53a6fbd9550071e..7ca671617832e6c6ee42577dfa48f673169887fc 100644 (file)
@@ -29,11 +29,17 @@ class AliMUONTrack : public TObject {
   AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegment, AliMUONEventReconstructor* EventReconstructor); // Constructor from two Segment's
   AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec, AliMUONEventReconstructor* EventReconstructor); // Constructor from one Segment and one HitForRec
 
+  AliMUONEventReconstructor* GetEventReconstructor(void) {return fEventReconstructor;};
   AliMUONTrackParam* GetTrackParamAtVertex(void);
-  void SetTrackParamAtVertex(void);
-  AliMUONTrackParam* GetTrackParamAtFirstHit(void);
+  void SetTrackParamAtVertex(void); // Set track parameters at vertex from last stations 4 & 5
+  void SetTrackParamAtVertex(AliMUONTrackParam* TrackParam) {fTrackParamAtVertex = *TrackParam;};
+
   TClonesArray* GetTrackHitsPtr(void);
   Int_t GetNTrackHits(void);
+  Int_t GetFitMCS(void) {return fFitMCS;};
+  void SetFitMCS(Int_t FitMCS);
+
+  AliMUONTrackParam* GetTrackParamAtFirstHit(void);
 
   void RecursiveDump(void); // Recursive dump (with track hits)
   void Fit(AliMUONTrackParam *TrackParam, Int_t NParam); // Fit
@@ -47,6 +53,7 @@ class AliMUONTrack : public TObject {
   AliMUONTrackParam fTrackParamAtVertex; // Track parameters at vertex
   TClonesArray *fTrackHitsPtr; // Pointer to array of TrackHit's
   Int_t fNTrackHits; // Number of TrackHit's
+  Int_t fFitMCS; // 0(1) for fit without(with) multiple Coulomb scattering
   
   ClassDef(AliMUONTrack, 1) // Class definition in ROOT context
     };
index e41ce1724fb7049289920f5fdc4f96219ffbeb36..1570db69d2a94571a3848bdbe4307432cecafcca 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.2  2000/06/15 07:58:49  morsch
+Code from MUON-dev joined
+
 Revision 1.1.2.3  2000/06/09 21:03:09  morsch
 Make includes consistent with new file structure.
 
@@ -229,3 +232,147 @@ void AliMUONTrackParam::ExtrapToStation(Int_t Station, AliMUONTrackParam *TrackP
   return;
 }
 
+  //__________________________________________________________________________
+void AliMUONTrackParam::ExtrapToVertex()
+{
+  // Extrapolation to the vertex.
+  // Returns the track parameters resulting from the extrapolation,
+  // in the current TrackParam.
+  // Changes parameters according to branson correction through the absorber 
+  
+  Double_t zAbsorber = 503.0; // to be coherent with the Geant absorber geometry !!!!
+  // Extrapolates track parameters upstream to the "Z" end of the front absorber
+  ExtrapToZ(zAbsorber);
+    // Makes Branson correction (multiple scattering + energy loss)
+  BransonCorrection();
+}
+
+  //__________________________________________________________________________
+void AliMUONTrackParam::BransonCorrection()
+{
+  // Branson correction of track parameters
+  // the entry parameters have to be calculated at the end of the absorber
+  Double_t zEndAbsorber, zBP, xBP, yBP;
+  Double_t  pYZ, pX, pY, pZ, pTotal, xEndAbsorber, yEndAbsorber, radiusEndAbsorber2, pT, theta;
+  Int_t sign;
+  // Would it be possible to calculate all that from Geant configuration ????
+  // Radiation lengths outer part theta > 3 degres
+  static Double_t x01[9] = { 18.8,    // C (cm)
+                            10.397,   // Concrete (cm)
+                            0.56,    // Plomb (cm)
+                            47.26,   // Polyethylene (cm)
+                            0.56,   // Plomb (cm)
+                            47.26,   // Polyethylene (cm)
+                            0.56,   // Plomb (cm)
+                            47.26,   // Polyethylene (cm)
+                            0.56 };   // Plomb (cm)
+  // inner part theta < 3 degres
+  static Double_t x02[3] = { 18.8,    // C (cm)
+                            10.397,   // Concrete (cm)
+                            0.35 };    // W (cm) 
+  // z positions of the materials inside the absober outer part theta > 3 degres
+  static Double_t z1[10] = { 90, 315, 467, 472, 477, 482, 487, 492, 497, 502 };
+  // inner part theta < 3 degres
+  static Double_t z2[4] = { 90, 315, 467, 503 };
+  static Bool_t first = kTRUE;
+  static Double_t zBP1, zBP2, rLimit;
+  // Calculates z positions of the Branson's planes: zBP1 for outer part and zBP2 for inner part (only at the first call)
+  if (first) {
+    first = kFALSE;
+    Double_t aNBP = 0.0;
+    Double_t aDBP = 0.0;
+    for (Int_t iBound = 0; iBound < 9; iBound++) {
+      aNBP = aNBP +
+       (z1[iBound+1] * z1[iBound+1] * z1[iBound+1] -
+        z1[iBound]   * z1[iBound]   * z1[iBound]    ) / x01[iBound];
+      aDBP = aDBP +
+       (z1[iBound+1] * z1[iBound+1] - z1[iBound]   * z1[iBound]    ) / x01[iBound];
+    }
+    zBP1 = (2.0 * aNBP) / (3.0 * aDBP);
+    aNBP = 0.0;
+    aDBP = 0.0;
+    for (Int_t iBound = 0; iBound < 3; iBound++) {
+      aNBP = aNBP +
+       (z2[iBound+1] * z2[iBound+1] * z2[iBound+1] -
+        z2[iBound]   * z2[iBound ]  * z2[iBound]    ) / x02[iBound];
+      aDBP = aDBP +
+       (z2[iBound+1] * z2[iBound+1] - z2[iBound] * z2[iBound]) / x02[iBound];
+    }
+    zBP2 = (2.0 * aNBP) / (3.0 * aDBP);
+    rLimit = z2[3] * TMath::Tan(3.0 * (TMath::Pi()) / 180.);
+  }
+
+  pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
+  sign = 1;      
+  if (fInverseBendingMomentum < 0) sign = -1;     
+  pZ = pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope)); 
+  pX = pZ * fNonBendingSlope; 
+  pY = pZ * fBendingSlope; 
+  pTotal = TMath::Sqrt(pYZ *pYZ + pX * pX);
+  xEndAbsorber = fNonBendingCoor; 
+  yEndAbsorber = fBendingCoor; 
+  radiusEndAbsorber2 = xEndAbsorber * xEndAbsorber + yEndAbsorber * yEndAbsorber;
+
+  if (radiusEndAbsorber2 > rLimit*rLimit) {
+    zEndAbsorber = z1[9];
+    zBP = zBP1;
+  } else {
+    zEndAbsorber = z2[2];
+    zBP = zBP2;
+  }
+
+  xBP = xEndAbsorber - (pX / pZ) * (zEndAbsorber - zBP);
+  yBP = yEndAbsorber - (pY / pZ) * (zEndAbsorber - zBP);
+
+  // new parameters after Branson and energy loss corrections
+  pZ = pTotal * zBP / TMath::Sqrt(xBP * xBP + yBP * yBP + zBP * zBP);
+  pX = pZ * xBP / zBP;
+  pY = pZ * yBP / zBP;
+  fBendingSlope = pY / pZ;
+  fNonBendingSlope = pX / pZ;
+  
+  pT = TMath::Sqrt(pX * pX + pY * pY);      
+  theta = TMath::ATan2(pT, pZ); 
+  pTotal =
+    TotalMomentumEnergyLoss(rLimit, pTotal, theta, xEndAbsorber, yEndAbsorber);
+
+  fInverseBendingMomentum = (sign / pTotal) *
+    TMath::Sqrt(1.0 +
+               fBendingSlope * fBendingSlope +
+               fNonBendingSlope * fNonBendingSlope) /
+    TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope);
+
+  // vertex position at (0,0,0)
+  // should be taken from vertex measurement ???
+  fBendingCoor = 0.0;
+  fNonBendingCoor = 0;
+  fZ= 0;
+}
+
+  //__________________________________________________________________________
+Double_t AliMUONTrackParam::TotalMomentumEnergyLoss(Double_t rLimit, Double_t pTotal, Double_t theta, Double_t xEndAbsorber, Double_t yEndAbsorber)
+{
+  // Returns the total momentum corrected from energy loss in the front absorber
+  Double_t deltaP, pTotalCorrected;
+
+  Double_t radiusEndAbsorber2 =
+    xEndAbsorber *xEndAbsorber + yEndAbsorber * yEndAbsorber;
+  // Parametrization to be redone according to change of absorber material ????
+  // The name is not so good, and there are many arguments !!!!
+  if (radiusEndAbsorber2 < rLimit * rLimit) {
+    if (pTotal < 15) {
+      deltaP = 2.737 + 0.0494 * pTotal - 0.001123 * pTotal * pTotal;
+    } else {
+      deltaP = 3.0643 + 0.01346 *pTotal;
+    }
+  } else {
+    if (pTotal < 15) {
+      deltaP  = 2.1380 + 0.0351 * pTotal - 0.000853 * pTotal * pTotal;
+    } else { 
+      deltaP = 2.407 + 0.00702 * pTotal;
+    }
+  }
+  pTotalCorrected = pTotal + deltaP / TMath::Cos(theta);
+  return pTotalCorrected;
+}
+
index 3ed93358fa9bee674276d916f39a1ca27d7a4f52..d089970d0b695f52290b541aaac8298687961278 100644 (file)
@@ -35,6 +35,9 @@ class AliMUONTrackParam : public TObject {
 
   void ExtrapToZ(Double_t Z);
   void ExtrapToStation(Int_t Station, AliMUONTrackParam *TrackParam);
+  void ExtrapToVertex();  // extrapolation to vertex through the absorber
+  void BransonCorrection(); // makes Branson correction
+  Double_t TotalMomentumEnergyLoss(Double_t rLimit, Double_t pTotal, Double_t theta, Double_t xEndAbsorber, Double_t yEndAbsorber); // returns total momentum after energy loss correction in the absorber
 
  protected:
  private: