[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$ */

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