[u/mrichter/AliRoot.git] / MUON / AliMUONTrackReconstructorK.cxx

/**************************************************************************
 * 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.                  *
 **************************************************************************/

/* $Id$ */

////////////////////////////////////
//
// MUON track reconstructor using the kalman method
//
// This class contains as data:
// * the parameters for the track reconstruction
//
// It contains as methods, among others:
// * MakeTracks to build the tracks
//
////////////////////////////////////

#include <stdlib.h>
#include <Riostream.h>
#include <TDirectory.h>
#include <TFile.h>
#include <TMatrixD.h>

#include "AliMUONVTrackReconstructor.h"
#include "AliMUONTrackReconstructorK.h"
#include "AliMUONData.h"
#include "AliMUONConstants.h"
#include "AliMUONHitForRec.h"
#include "AliMUONObjectPair.h"
#include "AliMUONRawCluster.h"
#include "AliMUONTrackK.h" 
#include "AliLog.h"

ClassImp(AliMUONTrackReconstructorK) // Class implementation in ROOT context

//__________________________________________________________________________
AliMUONTrackReconstructorK::AliMUONTrackReconstructorK(AliMUONData* data, const Option_t* TrackMethod)
  : AliMUONVTrackReconstructor(data),
    fTrackMethod(2), //tracking method (2-Kalman 3-Combination-Kalman/Clustering)
    fMuons(0)
{
  /// Constructor for class AliMUONTrackReconstructorK

  if (strstr(TrackMethod,"Kalman")) {
    cout << " *** Tracking with the Kalman filter *** " << endl;
    fTrackMethod = 2;
  } else if (strstr(TrackMethod,"Combi")) {
    cout << " *** Combined cluster / track finder ***" << endl;
    fTrackMethod = 3;
  } else AliFatal(Form("Tracking method %s not available",TrackMethod));
  
  // Memory allocation for the TClonesArray of reconstructed tracks
  fRecTracksPtr = new TClonesArray("AliMUONTrackK", 10);
}

  //__________________________________________________________________________
AliMUONTrackReconstructorK::~AliMUONTrackReconstructorK(void)
{
  /// Destructor for class AliMUONTrackReconstructorK
  delete fRecTracksPtr;
}

  //__________________________________________________________________________
void AliMUONTrackReconstructorK::AddHitsForRecFromRawClusters()
{
  /// To add to the list of hits for reconstruction all the raw clusters
  TTree *treeR;
  AliMUONHitForRec *hitForRec;
  AliMUONRawCluster *clus;
  Int_t iclus, nclus;
  TClonesArray *rawclusters;
  AliDebug(1,"Enter AddHitsForRecFromRawClusters");

  treeR = fMUONData->TreeR();
  if (!treeR) {
    AliError("TreeR must be loaded");
    exit(0);
  }
  
  if (fTrackMethod != 3) { //AZ
    fMUONData->SetTreeAddress("RC"); //AZ
    fMUONData->GetRawClusters(); // only one entry  
  }

  // Loop over tracking chambers
  for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
    rawclusters =fMUONData->RawClusters(ch);
    nclus = (Int_t) (rawclusters->GetEntries());
    // Loop over (cathode correlated) raw clusters
    for (iclus = 0; iclus < nclus; iclus++) {
      clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(iclus);
      // new AliMUONHitForRec from raw cluster
      // and increment number of AliMUONHitForRec's (total and in chamber)
      hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus);
      fNHitsForRec++;
      // more information into HitForRec
      hitForRec->SetBendingReso2(clus->GetErrY() * clus->GetErrY());
      hitForRec->SetNonBendingReso2(clus->GetErrX() * clus->GetErrX());
      //  original raw cluster
      hitForRec->SetChamberNumber(ch);
      hitForRec->SetHitNumber(iclus);
      // Z coordinate of the raw cluster (cm)
      hitForRec->SetZ(clus->GetZ(0));
      StdoutToAliDebug(3,
                       cout << "Chamber " << ch <<
                       " raw cluster  " << iclus << " : " << endl;
                       clus->Print("full");
                       cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
                       hitForRec->Print("full");
                       );
    } // end of cluster loop
  } // end of chamber loop
  SortHitsForRecWithIncreasingChamber(); 
  
  AliDebug(1,"End of AddHitsForRecFromRawClusters");
    if (AliLog::GetGlobalDebugLevel() > 0) {
      AliDebug(1, Form("NHitsForRec: %d",fNHitsForRec));
      for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
	AliDebug(1, Form("Chamber(0...): %d",ch));
	AliDebug(1, Form("NHitsForRec: %d", fNHitsForRecPerChamber[ch]));
	AliDebug(1, Form("Index(first HitForRec): %d", fIndexOfFirstHitForRecPerChamber[ch]));
	for (Int_t hit = fIndexOfFirstHitForRecPerChamber[ch];
	     hit < fIndexOfFirstHitForRecPerChamber[ch] + fNHitsForRecPerChamber[ch];
	     hit++) {
	  AliDebug(1, Form("HitForRec index(0...): %d",hit));
	  ((*fHitsForRecPtr)[hit])->Dump();
      }
    }
  }
  
  return;
}

  //__________________________________________________________________________
void AliMUONTrackReconstructorK::MakeTracks(void)
{
  /// To make the tracks from the list of segments and points in all stations
  AliDebug(1,"Enter MakeTracks");
  // The order may be important for the following Reset's
  //AZ ResetTracks();
  MakeTrackCandidates();
  if (fRecTracksPtr->GetEntriesFast() == 0) return;
  // Follow tracks in stations(1..) 3, 2 and 1
  FollowTracks();
  // Remove double tracks
  RemoveDoubleTracks();
  // Propagate tracks to the vertex through absorber
  ExtrapTracksToVertex();
  // Fill AliMUONTrack data members
  FillMUONTrack();
}

  //__________________________________________________________________________
void AliMUONTrackReconstructorK::MakeTrackCandidates(void)
{
  /// To make initial tracks for Kalman filter from segments in stations(1..)  4 and 5
  Int_t istat, iseg;
  TClonesArray *segments;
  AliMUONObjectPair *segment;
  AliMUONTrackK *trackK;

  AliDebug(1,"Enter MakeTrackCandidatesK");

  AliMUONTrackK a(this, fHitsForRecPtr);
  // Loop over stations(1...) 5 and 4
  for (istat=4; istat>=3; istat--) {
    // Make segments in the station
    segments = MakeSegmentsInStation(istat);
    // Loop over segments in the station
    for (iseg=0; iseg<segments->GetEntriesFast(); iseg++) {
      // Transform segments to tracks
      segment = (AliMUONObjectPair*) ((*segments)[iseg]);
      trackK = new ((*fRecTracksPtr)[fNRecTracks++]) AliMUONTrackK(segment);
    } // for (iseg=0;...)
  } // for (istat=4;...)
  return;
}

//__________________________________________________________________________
void AliMUONTrackReconstructorK::FollowTracks(void)
{
  /// Follow tracks using Kalman filter
  Bool_t ok;
  Int_t icand, ichamBeg = 0, ichamEnd, chamBits;
  AliMUONTrackK *trackK;
  AliMUONHitForRec *hit;
  AliMUONRawCluster *clus;
  TClonesArray *rawclusters;
  clus = 0; rawclusters = 0;

  Double_t simpleBPosition = 0.5 * (AliMUONConstants::CoilZ() + AliMUONConstants::YokeZ());
  Double_t simpleBLength = 0.5 * (AliMUONConstants::CoilL() + AliMUONConstants::YokeL());
  Double_t zDipole1 = simpleBPosition + 0.5 * simpleBLength;
  Double_t zDipole2 = zDipole1 - simpleBLength;

  // Print hits
  trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[0]);

  if (trackK->DebugLevel() > 0) {
    for (Int_t i1=0; i1<fNHitsForRec; i1++) {
      hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[i1]);
      printf(" Hit # %d %10.4f %10.4f %10.4f",
             hit->GetChamberNumber(), hit->GetBendingCoor(),
             hit->GetNonBendingCoor(), hit->GetZ());
 
      // from raw clusters
      rawclusters = fMUONData->RawClusters(hit->GetChamberNumber());
      clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->
							   GetHitNumber());
      printf(" %d", clus->GetTrack(1));
      if (clus->GetTrack(2) != -1) printf(" %d \n", clus->GetTrack(2));
      else printf("\n");
     
    }
  } // if (trackK->DebugLevel() > 0)

  icand = -1;
  Int_t nSeeds;
  nSeeds = fNRecTracks; // starting number of seeds
  // Loop over track candidates
  while (icand < fNRecTracks-1) {
    icand ++;
    if (trackK->DebugLevel()>0) cout << " *** Kalman track candidate No. " << icand << endl;
    trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
    if (trackK->GetRecover() < 0) continue; // failed track

    // Discard candidate which will produce the double track
    /*
    if (icand > 0) {
      ok = CheckCandidate(icand,nSeeds);
      if (!ok) {
        trackK->SetRecover(-1); // mark candidate to be removed
        continue;
      }
    }
    */

    ok = kTRUE;
    if (trackK->GetRecover() == 0) 
      hit = (AliMUONHitForRec*) trackK->GetTrackHits()->Last(); // last hit
    else 
      hit = trackK->GetHitLastOk(); // hit where track stopped

    if (hit) ichamBeg = hit->GetChamberNumber();
    ichamEnd = 0;
    // Check propagation direction
    if (!hit) { ichamBeg = ichamEnd; AliFatal(" ??? "); }
    else if (trackK->GetTrackDir() < 0) {
      ichamEnd = 9; // forward propagation
      ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
      if (ok) {
        ichamBeg = ichamEnd;
        ichamEnd = 6; // backward propagation
	// Change weight matrix and zero fChi2 for backpropagation
        trackK->StartBack();
	trackK->SetTrackDir(1);
        ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
        ichamBeg = ichamEnd;
        ichamEnd = 0;
      }
    } else {
      if (trackK->GetBPFlag()) {
	// backpropagation
        ichamEnd = 6; // backward propagation
	// Change weight matrix and zero fChi2 for backpropagation
        trackK->StartBack();
        ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
        ichamBeg = ichamEnd;
        ichamEnd = 0;
      }
    }

    if (ok) {
      trackK->SetTrackDir(1);
      trackK->SetBPFlag(kFALSE);
      ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
    }
    if (!ok) { trackK->SetRecover(-1); continue; } // mark candidate to be removed

    // Apply smoother
    if (trackK->GetRecover() >= 0) {
      ok = trackK->Smooth();
      if (!ok) trackK->SetRecover(-1); // mark candidate to be removed
    }

    // Majority 3 of 4 in first 2 stations
    if (!ok) continue;
    chamBits = 0;
    Double_t chi2max = 0;
    for (Int_t i=0; i<trackK->GetNTrackHits(); i++) {
      hit = (AliMUONHitForRec*) (*trackK->GetTrackHits())[i];
      chamBits |= BIT(hit->GetChamberNumber());
      if (trackK->GetChi2PerPoint(i) > chi2max) chi2max = trackK->GetChi2PerPoint(i);
    }
    if (!((chamBits&3)==3 || (chamBits>>2&3)==3) && chi2max > 25) {
      //trackK->Recover();
      trackK->SetRecover(-1); //mark candidate to be removed
      continue;
    }
    if (ok) trackK->SetTrackQuality(0); // compute "track quality"
  } // while

  for (Int_t i=0; i<fNRecTracks; i++) {
    trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
    if (trackK->GetRecover() < 0) fRecTracksPtr->RemoveAt(i);
  }

  // Compress TClonesArray
  fRecTracksPtr->Compress();
  fNRecTracks = fRecTracksPtr->GetEntriesFast();
  return;
}

//__________________________________________________________________________
Bool_t AliMUONTrackReconstructorK::CheckCandidate(Int_t icand, Int_t nSeeds) const
{
  /// Discards track candidate if it will produce the double track (having
  /// the same seed segment hits as hits of a good track found before)
  AliMUONTrackK *track1, *track2;
  AliMUONHitForRec *hit1, *hit2, *hit;
  AliMUONObjectPair *segment1, *segment2;

  track1 = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
  hit1 = (AliMUONHitForRec*) (*track1->GetTrackHits())[0]; // 1'st hit
  hit2 = (AliMUONHitForRec*) (*track1->GetTrackHits())[1]; // 2'nd hit

  for (Int_t i=0; i<icand; i++) {
    track2 = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
    //if (track2->GetRecover() < 0) continue;
    if (track2->GetRecover() < 0 && icand >= nSeeds) continue;
    segment1 = track1->GetStartSegment();
    segment2 = track2->GetStartSegment();
    if (segment1->First()  == segment2->First() &&
        segment1->Second() == segment2->Second()) {
      return kFALSE;
    } else {
      Int_t nSame = 0;
      for (Int_t j=0; j<track2->GetNTrackHits(); j++) {
        hit = (AliMUONHitForRec*) (*track2->GetTrackHits())[j];
        if (hit == hit1 || hit == hit2) {
          nSame++;
          if (nSame == 2) return kFALSE;
        }
      } // for (Int_t j=0;
    }
  } // for (Int_t i=0;
  return kTRUE;
}

  //__________________________________________________________________________
void AliMUONTrackReconstructorK::RemoveDoubleTracks(void)
{
  /// Removes double tracks (sharing more than half of their hits).
  /// Keeps the track with higher quality
  AliMUONTrackK *track1, *track2, *trackToKill;

  // Sort tracks according to their quality
  fRecTracksPtr->Sort();

  // Loop over first track of the pair
  track1 = (AliMUONTrackK*) fRecTracksPtr->First();
  Int_t debug = track1->DebugLevel();
  while (track1) {
    // Loop over second track of the pair
    track2 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
    while (track2) {
      // Check whether or not to keep track2
      if (!track2->KeepTrack(track1)) {
        if (debug >= 0) cout << " Killed track: " << 1/(*track2->GetTrackParameters())(4,0) <<
	  " " << track2->GetTrackQuality() << endl;
        trackToKill = track2;
        track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
        trackToKill->Kill();
        fRecTracksPtr->Compress();
      } else track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
    } // track2
    track1 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
  } // track1

  fNRecTracks = fRecTracksPtr->GetEntriesFast();
  if (debug >= 0) cout << " Number of Kalman tracks: " << fNRecTracks << endl;
}

  //__________________________________________________________________________
void AliMUONTrackReconstructorK::ExtrapTracksToVertex(void)
{
  /// Propagates track to the vertex thru absorber
  /// (using Branson correction for now)
  Double_t zVertex;
  zVertex = 0;
  for (Int_t i=0; i<fNRecTracks; i++) {
    //((AliMUONTrackK*)(*fRecTracksPtr)[i])->Branson();
    ((AliMUONTrackK*)(*fRecTracksPtr)[i])->SetTrackQuality(1); // compute Chi2
    //((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToZ(zVertex); // w/out absorber
    ((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToVertex(1); // with absorber
  }
}

  //__________________________________________________________________________
void AliMUONTrackReconstructorK::FillMUONTrack()
{
  // Set track parameters at hits for Kalman track. Fill fTrackParamAtHit of AliMUONTrack's
  AliMUONTrackK *track;
  track = (AliMUONTrackK*) fRecTracksPtr->First();
  while (track) {
    track->FillMUONTrack();
    track = (AliMUONTrackK*) fRecTracksPtr->After(track);
  } 
  return;
}

  //__________________________________________________________________________
void AliMUONTrackReconstructorK::EventDump(void)
{
  /// Dump reconstructed event (track parameters at vertex and at first hit),
  /// and the particle parameters
  AliDebug(1,"****** enter EventDump ******");
  AliDebug(1, Form("Number of Reconstructed tracks : %d", fNRecTracks)); 
  
  AliWarning("AliMUONTrackReconstructorK::EventDump: not implemented"); 
  
  // informations about generated particles NO !!!!!!!!
  
//    for (Int_t iPart = 0; iPart < np; iPart++) {
//      p = gAlice->Particle(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;
}
Commit	Line	Data
8cc77df0	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
	18	////////////////////////////////////
	19	//
	20	// MUON track reconstructor using the kalman method
	21	//
	22	// This class contains as data:
	23	// * the parameters for the track reconstruction
	24	//
	25	// It contains as methods, among others:
	26	// * MakeTracks to build the tracks
	27	//
	28	////////////////////////////////////
	29
	30	#include <stdlib.h>
	31	#include <Riostream.h>
	32	#include <TDirectory.h>
	33	#include <TFile.h>
	34	#include <TMatrixD.h>
	35
	36	#include "AliMUONVTrackReconstructor.h"
	37	#include "AliMUONTrackReconstructorK.h"
	38	#include "AliMUONData.h"
	39	#include "AliMUONConstants.h"
	40	#include "AliMUONHitForRec.h"
208f139e	41	#include "AliMUONObjectPair.h"
8cc77df0	42	#include "AliMUONRawCluster.h"
8cc77df0	43	#include "AliMUONTrackK.h"
	44	#include "AliLog.h"
	45
	46	ClassImp(AliMUONTrackReconstructorK) // Class implementation in ROOT context
	47
	48	//__________________________________________________________________________
	49	AliMUONTrackReconstructorK::AliMUONTrackReconstructorK(AliMUONData* data, const Option_t* TrackMethod)
	50	: AliMUONVTrackReconstructor(data),
	51	fTrackMethod(2), //tracking method (2-Kalman 3-Combination-Kalman/Clustering)
	52	fMuons(0)
	53	{
	54	/// Constructor for class AliMUONTrackReconstructorK
	55
	56	if (strstr(TrackMethod,"Kalman")) {
	57	cout << " * Tracking with the Kalman filter * " << endl;
	58	fTrackMethod = 2;
	59	} else if (strstr(TrackMethod,"Combi")) {
	60	cout << " * Combined cluster / track finder *" << endl;
	61	fTrackMethod = 3;
	62	} else AliFatal(Form("Tracking method %s not available",TrackMethod));
	63
	64	// Memory allocation for the TClonesArray of reconstructed tracks
	65	fRecTracksPtr = new TClonesArray("AliMUONTrackK", 10);
	66	}
	67
	68	//__________________________________________________________________________
	69	AliMUONTrackReconstructorK::~AliMUONTrackReconstructorK(void)
	70	{
	71	/// Destructor for class AliMUONTrackReconstructorK
	72	delete fRecTracksPtr;
	73	}
	74
	75	//__________________________________________________________________________
	76	void AliMUONTrackReconstructorK::AddHitsForRecFromRawClusters()
	77	{
	78	/// To add to the list of hits for reconstruction all the raw clusters
e6b25a6e	79	TTree *treeR;
8cc77df0	80	AliMUONHitForRec *hitForRec;
8cc77df0	81	AliMUONRawCluster *clus;
208f139e	82	Int_t iclus, nclus;
8cc77df0	83	TClonesArray *rawclusters;
	84	AliDebug(1,"Enter AddHitsForRecFromRawClusters");
	85
e6b25a6e	86	treeR = fMUONData->TreeR();
	87	if (!treeR) {
	88	AliError("TreeR must be loaded");
	89	exit(0);
	90	}
	91
8cc77df0	92	if (fTrackMethod != 3) { //AZ
e6b25a6e	93	fMUONData->SetTreeAddress("RC"); //AZ
8cc77df0	94	fMUONData->GetRawClusters(); // only one entry
	95	}
	96
	97	// Loop over tracking chambers
	98	for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
8cc77df0	99	rawclusters =fMUONData->RawClusters(ch);
	100	nclus = (Int_t) (rawclusters->GetEntries());
	101	// Loop over (cathode correlated) raw clusters
	102	for (iclus = 0; iclus < nclus; iclus++) {
	103	clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(iclus);
	104	// new AliMUONHitForRec from raw cluster
	105	// and increment number of AliMUONHitForRec's (total and in chamber)
	106	hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus);
	107	fNHitsForRec++;
8cc77df0	108	// more information into HitForRec
	109	hitForRec->SetBendingReso2(clus->GetErrY() * clus->GetErrY());
	110	hitForRec->SetNonBendingReso2(clus->GetErrX() * clus->GetErrX());
	111	// original raw cluster
	112	hitForRec->SetChamberNumber(ch);
	113	hitForRec->SetHitNumber(iclus);
	114	// Z coordinate of the raw cluster (cm)
	115	hitForRec->SetZ(clus->GetZ(0));
8cc77df0	116	StdoutToAliDebug(3,
	117	cout << "Chamber " << ch <<
	118	" raw cluster " << iclus << " : " << endl;
	119	clus->Print("full");
	120	cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
	121	hitForRec->Print("full");
	122	);
	123	} // end of cluster loop
	124	} // end of chamber loop
	125	SortHitsForRecWithIncreasingChamber();
	126
	127	AliDebug(1,"End of AddHitsForRecFromRawClusters");
	128	if (AliLog::GetGlobalDebugLevel() > 0) {
	129	AliDebug(1, Form("NHitsForRec: %d",fNHitsForRec));
	130	for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
	131	AliDebug(1, Form("Chamber(0...): %d",ch));
	132	AliDebug(1, Form("NHitsForRec: %d", fNHitsForRecPerChamber[ch]));
	133	AliDebug(1, Form("Index(first HitForRec): %d", fIndexOfFirstHitForRecPerChamber[ch]));
	134	for (Int_t hit = fIndexOfFirstHitForRecPerChamber[ch];
	135	hit < fIndexOfFirstHitForRecPerChamber[ch] + fNHitsForRecPerChamber[ch];
	136	hit++) {
	137	AliDebug(1, Form("HitForRec index(0...): %d",hit));
	138	((*fHitsForRecPtr)[hit])->Dump();
	139	}
	140	}
	141	}
	142
	143	return;
	144	}
	145
8cc77df0	146	//__________________________________________________________________________
	147	void AliMUONTrackReconstructorK::MakeTracks(void)
	148	{
	149	/// To make the tracks from the list of segments and points in all stations
	150	AliDebug(1,"Enter MakeTracks");
	151	// The order may be important for the following Reset's
	152	//AZ ResetTracks();
	153	MakeTrackCandidates();
	154	if (fRecTracksPtr->GetEntriesFast() == 0) return;
	155	// Follow tracks in stations(1..) 3, 2 and 1
	156	FollowTracks();
	157	// Remove double tracks
	158	RemoveDoubleTracks();
208f139e	159	// Propagate tracks to the vertex through absorber
208f139e	160	ExtrapTracksToVertex();
8cc77df0	161	// Fill AliMUONTrack data members
	162	FillMUONTrack();
	163	}
	164
	165	//__________________________________________________________________________
	166	void AliMUONTrackReconstructorK::MakeTrackCandidates(void)
	167	{
208f139e	168	/// To make initial tracks for Kalman filter from segments in stations(1..) 4 and 5
8cc77df0	169	Int_t istat, iseg;
208f139e	170	TClonesArray *segments;
208f139e	171	AliMUONObjectPair *segment;
8cc77df0	172	AliMUONTrackK *trackK;
	173
	174	AliDebug(1,"Enter MakeTrackCandidatesK");
	175
	176	AliMUONTrackK a(this, fHitsForRecPtr);
	177	// Loop over stations(1...) 5 and 4
	178	for (istat=4; istat>=3; istat--) {
208f139e	179	// Make segments in the station
208f139e	180	segments = MakeSegmentsInStation(istat);
8cc77df0	181	// Loop over segments in the station
208f139e	182	for (iseg=0; iseg<segments->GetEntriesFast(); iseg++) {
	183	// Transform segments to tracks
	184	segment = (AliMUONObjectPair) ((segments)[iseg]);
8cc77df0	185	trackK = new ((*fRecTracksPtr)[fNRecTracks++]) AliMUONTrackK(segment);
	186	} // for (iseg=0;...)
	187	} // for (istat=4;...)
	188	return;
	189	}
	190
	191	//__________________________________________________________________________
	192	void AliMUONTrackReconstructorK::FollowTracks(void)
	193	{
	194	/// Follow tracks using Kalman filter
	195	Bool_t ok;
	196	Int_t icand, ichamBeg = 0, ichamEnd, chamBits;
8cc77df0	197	AliMUONTrackK *trackK;
	198	AliMUONHitForRec *hit;
	199	AliMUONRawCluster *clus;
	200	TClonesArray *rawclusters;
	201	clus = 0; rawclusters = 0;
	202
208f139e	203	Double_t simpleBPosition = 0.5 * (AliMUONConstants::CoilZ() + AliMUONConstants::YokeZ());
	204	Double_t simpleBLength = 0.5 * (AliMUONConstants::CoilL() + AliMUONConstants::YokeL());
	205	Double_t zDipole1 = simpleBPosition + 0.5 * simpleBLength;
	206	Double_t zDipole2 = zDipole1 - simpleBLength;
8cc77df0	207
	208	// Print hits
	209	trackK = (AliMUONTrackK) ((fRecTracksPtr)[0]);
	210
	211	if (trackK->DebugLevel() > 0) {
	212	for (Int_t i1=0; i1<fNHitsForRec; i1++) {
	213	hit = (AliMUONHitForRec) ((fHitsForRecPtr)[i1]);
	214	printf(" Hit # %d %10.4f %10.4f %10.4f",
	215	hit->GetChamberNumber(), hit->GetBendingCoor(),
	216	hit->GetNonBendingCoor(), hit->GetZ());
	217
	218	// from raw clusters
	219	rawclusters = fMUONData->RawClusters(hit->GetChamberNumber());
	220	clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->
	221	GetHitNumber());
	222	printf(" %d", clus->GetTrack(1));
	223	if (clus->GetTrack(2) != -1) printf(" %d \n", clus->GetTrack(2));
	224	else printf("\n");
	225
	226	}
	227	} // if (trackK->DebugLevel() > 0)
	228
	229	icand = -1;
	230	Int_t nSeeds;
	231	nSeeds = fNRecTracks; // starting number of seeds
	232	// Loop over track candidates
	233	while (icand < fNRecTracks-1) {
	234	icand ++;
	235	if (trackK->DebugLevel()>0) cout << " *** Kalman track candidate No. " << icand << endl;
	236	trackK = (AliMUONTrackK) ((fRecTracksPtr)[icand]);
	237	if (trackK->GetRecover() < 0) continue; // failed track
	238
	239	// Discard candidate which will produce the double track
	240	/*
	241	if (icand > 0) {
	242	ok = CheckCandidate(icand,nSeeds);
	243	if (!ok) {
	244	trackK->SetRecover(-1); // mark candidate to be removed
	245	continue;
	246	}
	247	}
	248	*/
	249
	250	ok = kTRUE;
	251	if (trackK->GetRecover() == 0)
	252	hit = (AliMUONHitForRec*) trackK->GetTrackHits()->Last(); // last hit
	253	else
	254	hit = trackK->GetHitLastOk(); // hit where track stopped
	255
	256	if (hit) ichamBeg = hit->GetChamberNumber();
	257	ichamEnd = 0;
	258	// Check propagation direction
	259	if (!hit) { ichamBeg = ichamEnd; AliFatal(" ??? "); }
	260	else if (trackK->GetTrackDir() < 0) {
	261	ichamEnd = 9; // forward propagation
	262	ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
	263	if (ok) {
	264	ichamBeg = ichamEnd;
	265	ichamEnd = 6; // backward propagation
	266	// Change weight matrix and zero fChi2 for backpropagation
	267	trackK->StartBack();
	268	trackK->SetTrackDir(1);
	269	ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
	270	ichamBeg = ichamEnd;
271	ichamEnd = 0;
272	}
273	} else {
274	if (trackK->GetBPFlag()) {
275	// backpropagation
276	ichamEnd = 6; // backward propagation
277	// Change weight matrix and zero fChi2 for backpropagation
278	trackK->StartBack();
279	ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
280	ichamBeg = ichamEnd;
281	ichamEnd = 0;
282	}
283	}
284
285	if (ok) {
286	trackK->SetTrackDir(1);
287	trackK->SetBPFlag(kFALSE);
288	ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
289	}
290	if (!ok) { trackK->SetRecover(-1); continue; } // mark candidate to be removed
291
292	// Apply smoother
293	if (trackK->GetRecover() >= 0) {
294	ok = trackK->Smooth();
295	if (!ok) trackK->SetRecover(-1); // mark candidate to be removed
296	}
297
298	// Majority 3 of 4 in first 2 stations
299	if (!ok) continue;
300	chamBits = 0;
301	Double_t chi2max = 0;
302	for (Int_t i=0; i<trackK->GetNTrackHits(); i++) {
303	hit = (AliMUONHitForRec) (trackK->GetTrackHits())[i];
304	chamBits \|= BIT(hit->GetChamberNumber());
305	if (trackK->GetChi2PerPoint(i) > chi2max) chi2max = trackK->GetChi2PerPoint(i);
306	}
307	if (!((chamBits&3)==3 \|\| (chamBits>>2&3)==3) && chi2max > 25) {
308	//trackK->Recover();
309	trackK->SetRecover(-1); //mark candidate to be removed
310	continue;
311	}
312	if (ok) trackK->SetTrackQuality(0); // compute "track quality"
313	} // while
314
315	for (Int_t i=0; i<fNRecTracks; i++) {
316	trackK = (AliMUONTrackK) ((fRecTracksPtr)[i]);
317	if (trackK->GetRecover() < 0) fRecTracksPtr->RemoveAt(i);
318	}
319
320	// Compress TClonesArray
321	fRecTracksPtr->Compress();
322	fNRecTracks = fRecTracksPtr->GetEntriesFast();
323	return;
324	}
325
326	//__________________________________________________________________________
327	Bool_t AliMUONTrackReconstructorK::CheckCandidate(Int_t icand, Int_t nSeeds) const
328	{
329	/// Discards track candidate if it will produce the double track (having
330	/// the same seed segment hits as hits of a good track found before)
331	AliMUONTrackK track1, track2;
332	AliMUONHitForRec hit1, hit2, *hit;
208f139e	333	AliMUONObjectPair segment1, segment2;
8cc77df0	334
	335	track1 = (AliMUONTrackK) ((fRecTracksPtr)[icand]);
	336	hit1 = (AliMUONHitForRec) (track1->GetTrackHits())[0]; // 1'st hit
	337	hit2 = (AliMUONHitForRec) (track1->GetTrackHits())[1]; // 2'nd hit
	338
	339	for (Int_t i=0; i<icand; i++) {
	340	track2 = (AliMUONTrackK) ((fRecTracksPtr)[i]);
	341	//if (track2->GetRecover() < 0) continue;
	342	if (track2->GetRecover() < 0 && icand >= nSeeds) continue;
208f139e	343	segment1 = track1->GetStartSegment();
	344	segment2 = track2->GetStartSegment();
	345	if (segment1->First() == segment2->First() &&
	346	segment1->Second() == segment2->Second()) {
8cc77df0	347	return kFALSE;
	348	} else {
	349	Int_t nSame = 0;
	350	for (Int_t j=0; j<track2->GetNTrackHits(); j++) {
	351	hit = (AliMUONHitForRec) (track2->GetTrackHits())[j];
	352	if (hit == hit1 \|\| hit == hit2) {
	353	nSame++;
	354	if (nSame == 2) return kFALSE;
	355	}
	356	} // for (Int_t j=0;
	357	}
	358	} // for (Int_t i=0;
	359	return kTRUE;
	360	}
	361
	362	//__________________________________________________________________________
	363	void AliMUONTrackReconstructorK::RemoveDoubleTracks(void)
	364	{
	365	/// Removes double tracks (sharing more than half of their hits).
	366	/// Keeps the track with higher quality
	367	AliMUONTrackK track1, track2, *trackToKill;
	368
	369	// Sort tracks according to their quality
	370	fRecTracksPtr->Sort();
	371
	372	// Loop over first track of the pair
	373	track1 = (AliMUONTrackK*) fRecTracksPtr->First();
	374	Int_t debug = track1->DebugLevel();
	375	while (track1) {
	376	// Loop over second track of the pair
	377	track2 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
	378	while (track2) {
	379	// Check whether or not to keep track2
	380	if (!track2->KeepTrack(track1)) {
	381	if (debug >= 0) cout << " Killed track: " << 1/(*track2->GetTrackParameters())(4,0) <<
	382	" " << track2->GetTrackQuality() << endl;
	383	trackToKill = track2;
	384	track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
	385	trackToKill->Kill();
	386	fRecTracksPtr->Compress();
	387	} else track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
	388	} // track2
	389	track1 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
	390	} // track1
	391
	392	fNRecTracks = fRecTracksPtr->GetEntriesFast();
	393	if (debug >= 0) cout << " Number of Kalman tracks: " << fNRecTracks << endl;
	394	}
	395
	396	//__________________________________________________________________________
208f139e	397	void AliMUONTrackReconstructorK::ExtrapTracksToVertex(void)
8cc77df0	398	{
	399	/// Propagates track to the vertex thru absorber
	400	/// (using Branson correction for now)
	401	Double_t zVertex;
	402	zVertex = 0;
	403	for (Int_t i=0; i<fNRecTracks; i++) {
	404	//((AliMUONTrackK)(fRecTracksPtr)[i])->Branson();
	405	((AliMUONTrackK)(fRecTracksPtr)[i])->SetTrackQuality(1); // compute Chi2
	406	//((AliMUONTrackK)(fRecTracksPtr)[i])->GoToZ(zVertex); // w/out absorber
	407	((AliMUONTrackK)(fRecTracksPtr)[i])->GoToVertex(1); // with absorber
	408	}
	409	}
	410
	411	//__________________________________________________________________________
	412	void AliMUONTrackReconstructorK::FillMUONTrack()
	413	{
	414	// Set track parameters at hits for Kalman track. Fill fTrackParamAtHit of AliMUONTrack's
	415	AliMUONTrackK *track;
	416	track = (AliMUONTrackK*) fRecTracksPtr->First();
	417	while (track) {
	418	track->FillMUONTrack();
	419	track = (AliMUONTrackK*) fRecTracksPtr->After(track);
	420	}
	421	return;
	422	}
	423
	424	//__________________________________________________________________________
	425	void AliMUONTrackReconstructorK::EventDump(void)
	426	{
	427	/// Dump reconstructed event (track parameters at vertex and at first hit),
	428	/// and the particle parameters
	429	AliDebug(1,"**** enter EventDump ****");
	430	AliDebug(1, Form("Number of Reconstructed tracks : %d", fNRecTracks));
	431
	432	AliWarning("AliMUONTrackReconstructorK::EventDump: not implemented");
	433
	434	// informations about generated particles NO !!!!!!!!
	435
	436	// for (Int_t iPart = 0; iPart < np; iPart++) {
	437	// p = gAlice->Particle(iPart);
	438	// printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n",
	439	// iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode());
	440	// }
	441	return;
	442	}
	443
	444