[u/mrichter/AliRoot.git] / MUON / AliMUONRecoCheck.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 AliMUONRecoCheck
// -----------------------
// Utility class to check the muon reconstruction. Reconstructed tracks are compared
// to reference tracks. The reference tracks are built from AliTrackReference for the
// hit in chamber (0..9) and from kinematics for the vertex parameters.     

#include <TParticle.h>

#include "AliRun.h" // for gAlice
#include "AliLog.h" 
#include "AliLoader.h" 
#include "AliRunLoader.h" 
#include "AliTrackReference.h"
#include "AliHeader.h"
#include "AliMC.h"
#include "AliStack.h"
#include "AliMUON.h"
#include "AliMUONRecoCheck.h"
#include "AliMUONTrack.h"
#include "AliMUONData.h"
#include "AliMUONConstants.h"

ClassImp(AliMUONRecoCheck)

//_____________________________________________________________________________
AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader)
{
/// Constructor
  fMuonTrackRef = new TClonesArray("AliMUONTrack", 10);

  // open the run loader
  fRunLoader = AliRunLoader::Open(chLoader);
  if (!fRunLoader) {
    AliError(Form("no run loader found in file %s","galice.root" ));
    return;
  }
  // initialize loader's
  AliLoader *loader = fRunLoader->GetLoader("MUONLoader");

  // initialize container
  fMUONData  = new AliMUONData(loader,"MUON","MUON");

   // Loading AliRun master
  if (fRunLoader->GetAliRun() == 0x0) fRunLoader->LoadgAlice();

  fRunLoader->LoadKinematics("READ");
  fRunLoader->LoadTrackRefs("READ");
  loader->LoadTracks("READ");

  fReconstructibleTracks = 0; 
  fRecoTracks = 0;
}

//____________________________________________________________________
AliMUONRecoCheck::AliMUONRecoCheck(const AliMUONRecoCheck& rhs)
 : TObject(rhs)
{
/// Protected copy constructor

  AliFatal("Not implemented.");
}

//_____________________________________________________________________________
AliMUONRecoCheck::~AliMUONRecoCheck()
{
/// Destructor

  fRunLoader->UnloadKinematics();
  fRunLoader->UnloadTrackRefs();
  fRunLoader->UnloadTracks();
  fMuonTrackRef->Delete();
  delete fMuonTrackRef;
  delete fMUONData;
}

//________________________________________________________________________
AliMUONRecoCheck& AliMUONRecoCheck::operator = (const AliMUONRecoCheck& rhs)
{
/// Protected assignement operator

  if (this == &rhs) return *this;

  AliFatal("Not implemented.");
    
  return *this;  
}

//_____________________________________________________________________________
void AliMUONRecoCheck::MakeTrackRef()
{
/// Make reconstructible tracks

  AliTrackReference *trackReference;
  AliMUONTrack *muonTrack;
  AliMUONTrackParam *trackParam;
  AliMUONHitForRec *hitForRec;
  Float_t x, y, z, pX, pY, pZ, pYZ;
  Int_t track, trackSave;
  TParticle *particle;   
  Float_t bendingSlope = 0;
  Float_t  nonBendingSlope = 0;
  Float_t inverseBendingMomentum = 0;
  
  TTree* treeTR = fRunLoader->TreeTR();
  if (treeTR == NULL) return;

  TBranch* branch = treeTR->GetBranch("MUON");
  if (branch == NULL) return;

  TClonesArray* trackRefs = 0;
  branch->SetAddress(&trackRefs);
  branch->SetAutoDelete(kTRUE);  

  Int_t nTrackRef = (Int_t)branch->GetEntries();
 
  track = trackSave = -999;
  Bool_t isNewTrack;
  Int_t iHitMin, iChamber;

  trackParam = new AliMUONTrackParam();
  hitForRec = new AliMUONHitForRec();
  muonTrack = new AliMUONTrack();

  for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; iTrackRef++) {
    branch->GetEntry(iTrackRef);
    
    iHitMin = 0;
    isNewTrack = kTRUE;
    
    if (!trackRefs->GetEntries()) continue; 

    while (isNewTrack) {
      
      for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); iHit++) {
      
	trackReference = (AliTrackReference*)trackRefs->At(iHit);
	x = trackReference->X();
	y = trackReference->Y();
	z = trackReference->Z();
	pX = trackReference->Px();
	pY = trackReference->Py();
	pZ = trackReference->Pz();

	track = trackReference->GetTrack();

	if (track != trackSave && iHit != 0) {
	  iHitMin = iHit;
	  trackSave = track;
	  break;
	}

	// track parameters at hit
	trackParam->SetBendingCoor(y);
	trackParam->SetNonBendingCoor(x);
	trackParam->SetZ(z);
       
	if (TMath::Abs(pZ) > 0) {
	  bendingSlope = pY/pZ;
	  nonBendingSlope = pX/pZ;
	}
	pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
	if (pYZ >0) inverseBendingMomentum = 1/pYZ; 

	trackParam->SetBendingSlope(bendingSlope);
	trackParam->SetNonBendingSlope(nonBendingSlope);
	trackParam->SetInverseBendingMomentum(inverseBendingMomentum);

	hitForRec->SetBendingCoor(y);
	hitForRec->SetNonBendingCoor(x);
	hitForRec->SetZ(z);
	hitForRec->SetBendingReso2(0.0); 
	hitForRec->SetNonBendingReso2(0.0);  
	iChamber = AliMUONConstants::ChamberNumber(z);
	hitForRec->SetChamberNumber(iChamber);

	muonTrack->AddTrackParamAtHit(trackParam);
	muonTrack->AddHitForRecAtHit(hitForRec);
	muonTrack->SetTrackID(track);
	
	trackSave = track;
	if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
      }
      
      // track parameters at vertex 
      particle = fRunLoader->GetHeader()->Stack()->Particle(muonTrack->GetTrackID());
      if (particle) {
	x = particle->Vx();
	y = particle->Vy();
	z = particle->Vz();
	pX = particle->Px();
	pY = particle->Py();
	pZ = particle->Pz();

	trackParam->SetBendingCoor(y);
	trackParam->SetNonBendingCoor(x);
	trackParam->SetZ(z);
	if (TMath::Abs(pZ) > 0) {
	  bendingSlope = pY/pZ;
	  nonBendingSlope = pX/pZ;
	}
	pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
	if (pYZ >0) inverseBendingMomentum = 1/pYZ;       
	trackParam->SetBendingSlope(bendingSlope);
	trackParam->SetNonBendingSlope(nonBendingSlope);
	trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
      
	muonTrack->SetTrackParamAtVertex(trackParam);
      }
      
      AddMuonTrackReference(muonTrack);
      muonTrack->ResetTrackParamAtHit();
      muonTrack->ResetHitForRecAtHit();
      
    } // end while isNewTrack
    
  }
  
  CleanMuonTrackRef();
  
  ReconstructibleTracks();

  delete muonTrack;
  delete trackParam;
  delete hitForRec;

}

//____________________________________________________________________________
TClonesArray* AliMUONRecoCheck::GetTrackReco()
{
/// Return TClonesArray of reconstructed tracks

  GetMUONData()->ResetRecTracks();
  GetMUONData()->SetTreeAddress("RT");
  fTrackReco = GetMUONData()->RecTracks(); 
  GetMUONData()->GetRecTracks();
  fRecoTracks = fTrackReco->GetEntriesFast();
  return fTrackReco;
}
//_____________________________________________________________________________
void AliMUONRecoCheck::PrintEvent() const
{
/// Debug facility

  AliMUONTrack *track;
  AliMUONHitForRec *hitForRec;
  TClonesArray *  hitForRecAtHit = 0;
  Float_t xRec,yRec,zRec;

  Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
  
  printf(" ******************************* \n");
  printf(" nb of tracks %d \n",nTrackRef);

  for (Int_t index = 0; index < nTrackRef; index++) {
    track = (AliMUONTrack*)fMuonTrackRef->At(index);
    hitForRecAtHit = track->GetHitForRecAtHit();
    Int_t nTrackHits = hitForRecAtHit->GetEntriesFast();
    printf(" track number %d \n",index);
    for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
      hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit); 
      xRec  = hitForRec->GetNonBendingCoor();
      yRec  = hitForRec->GetBendingCoor();
      zRec  = hitForRec->GetZ();
      printf("  x,y,z: %f , %f , %f \n",xRec,yRec,zRec);
    }
  }
}

//_____________________________________________________________________________
void AliMUONRecoCheck::ResetTracks() const
{
/// Reset tracks

  if (fMuonTrackRef) fMuonTrackRef->Clear();
}
//_____________________________________________________________________________
void AliMUONRecoCheck::CleanMuonTrackRef()
{
/// Re-calculate hits parameters because two AliTrackReferences are recorded for
/// each chamber (one when particle is entering + one when particle is leaving 
/// the sensitive volume) 

  Float_t maxGasGap = 1.; // cm 
  AliMUONTrack *track, *trackNew;
  AliMUONHitForRec *hitForRec, *hitForRec1, *hitForRec2;
  AliMUONTrackParam *trackParam, *trackParam1, *trackParam2, *trackParamAtVertex;
  TClonesArray *  hitForRecAtHit = 0;
  TClonesArray *  trackParamAtHit = 0;
  Float_t xRec,yRec,zRec;
  Float_t xRec1,yRec1,zRec1;
  Float_t xRec2,yRec2,zRec2;
  Float_t bendingSlope,nonBendingSlope,bendingMomentum;
  Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
  Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
  TClonesArray *newMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
  Int_t iHit1;
  Int_t iChamber = 0;
  Int_t nRec = 0;
  Int_t nTrackHits = 0;

  hitForRec = new AliMUONHitForRec();
  trackParam = new AliMUONTrackParam();
  trackNew = new AliMUONTrack();

  Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
  
  for (Int_t index = 0; index < nTrackRef; index++) {
    track = (AliMUONTrack*)fMuonTrackRef->At(index);
    hitForRecAtHit = track->GetHitForRecAtHit();
    trackParamAtHit = track->GetTrackParamAtHit();
    trackParamAtVertex = track->GetTrackParamAtVertex();
    nTrackHits = hitForRecAtHit->GetEntriesFast();
    iHit1 = 0;
    while (iHit1 < nTrackHits) {
      hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1); 
      trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1); 
      xRec1  = hitForRec1->GetNonBendingCoor();
      yRec1  = hitForRec1->GetBendingCoor();
      zRec1  = hitForRec1->GetZ();	
      xRec   = xRec1;
      yRec   = yRec1;
      zRec   = zRec1;
      bendingSlope1 = trackParam1->GetBendingSlope();
      nonBendingSlope1 = trackParam1->GetNonBendingSlope();
      bendingMomentum1 = 0;
      if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
	bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
      bendingSlope = bendingSlope1;
      nonBendingSlope = nonBendingSlope1;
      bendingMomentum = bendingMomentum1;
      nRec = 1;  
      for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
	hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2); 
	trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2); 
	xRec2  = hitForRec2->GetNonBendingCoor();
	yRec2  = hitForRec2->GetBendingCoor();
	zRec2  = hitForRec2->GetZ();	  
	bendingSlope2 = trackParam2->GetBendingSlope();
	nonBendingSlope2 = trackParam2->GetNonBendingSlope();
	bendingMomentum2 = 0;
	if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
	  bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
	
	if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
	  nRec++;
	  xRec += xRec2;
	  yRec += yRec2;
	  zRec += zRec2;
	  bendingSlope += bendingSlope2;
	  nonBendingSlope += nonBendingSlope2;
	  bendingMomentum += bendingMomentum2;
	  iHit1 = iHit2;
	}
	
      } // end iHit2
      xRec /= (Float_t)nRec;
      yRec /= (Float_t)nRec;
      zRec /= (Float_t)nRec;
      bendingSlope /= (Float_t)nRec;
      nonBendingSlope /= (Float_t)nRec;
      bendingMomentum /= (Float_t)nRec;
      
      hitForRec->SetNonBendingCoor(xRec);
      hitForRec->SetBendingCoor(yRec);
      hitForRec->SetZ(zRec);
      iChamber = AliMUONConstants::ChamberNumber(zRec);
      hitForRec->SetChamberNumber(iChamber);
      hitForRec->SetBendingReso2(0.0); 
      hitForRec->SetNonBendingReso2(0.0); 
      trackParam->SetNonBendingCoor(xRec);
      trackParam->SetBendingCoor(yRec);
      trackParam->SetZ(zRec);
      trackParam->SetNonBendingSlope(nonBendingSlope);
      trackParam->SetBendingSlope(bendingSlope);
      if (TMath::Abs(bendingMomentum) > 0)
	trackParam->SetInverseBendingMomentum(1./bendingMomentum);

      trackNew->AddHitForRecAtHit(hitForRec);
      trackNew->AddTrackParamAtHit(trackParam);
      
      iHit1++;
    } // end iHit1

    trackNew->SetTrackID(track->GetTrackID());
    trackNew->SetTrackParamAtVertex(trackParamAtVertex);
    {new ((*newMuonTrackRef)[newMuonTrackRef->GetEntriesFast()]) AliMUONTrack(*trackNew);}    
    trackNew->ResetHitForRecAtHit();
    trackNew->ResetTrackParamAtHit();
    
  } // end trackRef

  fMuonTrackRef->Clear();
  nTrackRef = newMuonTrackRef->GetEntriesFast();
  for (Int_t index = 0; index < nTrackRef; index++) {
    track = (AliMUONTrack*)newMuonTrackRef->At(index);
    AddMuonTrackReference(track);
  }
  

  delete trackNew;
  delete hitForRec;
  delete trackParam;
  newMuonTrackRef->Delete();
  delete newMuonTrackRef;
  
}

//_____________________________________________________________________________
void AliMUONRecoCheck::ReconstructibleTracks()
{
/// Calculate the number of reconstructible tracks
  
  AliMUONTrack* track;
  TClonesArray* hitForRecAtHit = NULL;
  AliMUONHitForRec* hitForRec;
  Float_t zRec;
  Int_t nTrackRef, nTrackHits;
  Int_t isChamberInTrack[10];
  Int_t iChamber = 0;
  Bool_t isTrackOK = kTRUE;

  fReconstructibleTracks = 0;

  nTrackRef = fMuonTrackRef->GetEntriesFast();
  for (Int_t index = 0; index < nTrackRef; index++) {
    track = (AliMUONTrack*)fMuonTrackRef->At(index);
    hitForRecAtHit = track->GetHitForRecAtHit();
    nTrackHits = hitForRecAtHit->GetEntriesFast();
    for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
 
    for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
      hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit); 
      zRec  = hitForRec->GetZ();
      iChamber = hitForRec->GetChamberNumber();
      if (iChamber < 0 || iChamber > 10) continue;
      isChamberInTrack[iChamber] = 1;
      
    }
    // track is reconstructible if the particle is crossing every tracking chambers
    isTrackOK = kTRUE;
    for (Int_t ch = 0; ch < 10; ch++) {
      if (!isChamberInTrack[ch]) isTrackOK = kFALSE;
    }    if (isTrackOK) fReconstructibleTracks++;
    if (!isTrackOK) fMuonTrackRef->Remove(track); // remove non reconstructible tracks
  }
  fMuonTrackRef->Compress();
}
Commit	Line	Data
	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	// Class AliMUONRecoCheck
	20	// -----------------------
	21	// Utility class to check the muon reconstruction. Reconstructed tracks are compared
	22	// to reference tracks. The reference tracks are built from AliTrackReference for the
	23	// hit in chamber (0..9) and from kinematics for the vertex parameters.
	24
	25	#include <TParticle.h>
	26
	27	#include "AliRun.h" // for gAlice
	28	#include "AliLog.h"
	29	#include "AliLoader.h"
	30	#include "AliRunLoader.h"
	31	#include "AliTrackReference.h"
	32	#include "AliHeader.h"
	33	#include "AliMC.h"
	34	#include "AliStack.h"
	35	#include "AliMUON.h"
	36	#include "AliMUONRecoCheck.h"
	37	#include "AliMUONTrack.h"
	38	#include "AliMUONData.h"
	39	#include "AliMUONConstants.h"
	40
	41	ClassImp(AliMUONRecoCheck)
	42
	43	//_____________________________________________________________________________
	44	AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader)
	45	{
	46	/// Constructor
	47	fMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
	48
	49	// open the run loader
	50	fRunLoader = AliRunLoader::Open(chLoader);
	51	if (!fRunLoader) {
	52	AliError(Form("no run loader found in file %s","galice.root" ));
	53	return;
	54	}
	55	// initialize loader's
	56	AliLoader *loader = fRunLoader->GetLoader("MUONLoader");
	57
	58	// initialize container
	59	fMUONData = new AliMUONData(loader,"MUON","MUON");
	60
	61	// Loading AliRun master
	62	if (fRunLoader->GetAliRun() == 0x0) fRunLoader->LoadgAlice();
	63
	64	fRunLoader->LoadKinematics("READ");
	65	fRunLoader->LoadTrackRefs("READ");
	66	loader->LoadTracks("READ");
	67
	68	fReconstructibleTracks = 0;
	69	fRecoTracks = 0;
	70	}
	71
	72	//____________________________________________________________________
	73	AliMUONRecoCheck::AliMUONRecoCheck(const AliMUONRecoCheck& rhs)
	74	: TObject(rhs)
	75	{
	76	/// Protected copy constructor
	77
	78	AliFatal("Not implemented.");
	79	}
	80
	81	//_____________________________________________________________________________
	82	AliMUONRecoCheck::~AliMUONRecoCheck()
	83	{
	84	/// Destructor
	85
	86	fRunLoader->UnloadKinematics();
	87	fRunLoader->UnloadTrackRefs();
	88	fRunLoader->UnloadTracks();
	89	fMuonTrackRef->Delete();
	90	delete fMuonTrackRef;
	91	delete fMUONData;
	92	}
	93
	94	//________________________________________________________________________
	95	AliMUONRecoCheck& AliMUONRecoCheck::operator = (const AliMUONRecoCheck& rhs)
	96	{
	97	/// Protected assignement operator
	98
	99	if (this == &rhs) return *this;
	100
	101	AliFatal("Not implemented.");
	102
	103	return *this;
	104	}
	105
	106	//_____________________________________________________________________________
	107	void AliMUONRecoCheck::MakeTrackRef()
	108	{
	109	/// Make reconstructible tracks
	110
	111	AliTrackReference *trackReference;
	112	AliMUONTrack *muonTrack;
	113	AliMUONTrackParam *trackParam;
	114	AliMUONHitForRec *hitForRec;
	115	Float_t x, y, z, pX, pY, pZ, pYZ;
	116	Int_t track, trackSave;
	117	TParticle *particle;
	118	Float_t bendingSlope = 0;
	119	Float_t nonBendingSlope = 0;
	120	Float_t inverseBendingMomentum = 0;
	121
	122	TTree* treeTR = fRunLoader->TreeTR();
	123	if (treeTR == NULL) return;
	124
	125	TBranch* branch = treeTR->GetBranch("MUON");
	126	if (branch == NULL) return;
	127
	128	TClonesArray* trackRefs = 0;
	129	branch->SetAddress(&trackRefs);
	130	branch->SetAutoDelete(kTRUE);
	131
	132	Int_t nTrackRef = (Int_t)branch->GetEntries();
	133
	134	track = trackSave = -999;
	135	Bool_t isNewTrack;
	136	Int_t iHitMin, iChamber;
	137
	138	trackParam = new AliMUONTrackParam();
	139	hitForRec = new AliMUONHitForRec();
	140	muonTrack = new AliMUONTrack();
	141
	142	for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; iTrackRef++) {
	143	branch->GetEntry(iTrackRef);
	144
	145	iHitMin = 0;
	146	isNewTrack = kTRUE;
	147
	148	if (!trackRefs->GetEntries()) continue;
	149
	150	while (isNewTrack) {
	151
	152	for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); iHit++) {
	153
	154	trackReference = (AliTrackReference*)trackRefs->At(iHit);
	155	x = trackReference->X();
	156	y = trackReference->Y();
	157	z = trackReference->Z();
	158	pX = trackReference->Px();
	159	pY = trackReference->Py();
	160	pZ = trackReference->Pz();
	161
	162	track = trackReference->GetTrack();
	163
	164	if (track != trackSave && iHit != 0) {
	165	iHitMin = iHit;
	166	trackSave = track;
	167	break;
	168	}
	169
	170	// track parameters at hit
	171	trackParam->SetBendingCoor(y);
	172	trackParam->SetNonBendingCoor(x);
	173	trackParam->SetZ(z);
	174
	175	if (TMath::Abs(pZ) > 0) {
	176	bendingSlope = pY/pZ;
	177	nonBendingSlope = pX/pZ;
	178	}
	179	pYZ = TMath::Sqrt(pYpY+pZpZ);
	180	if (pYZ >0) inverseBendingMomentum = 1/pYZ;
	181
	182	trackParam->SetBendingSlope(bendingSlope);
	183	trackParam->SetNonBendingSlope(nonBendingSlope);
	184	trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
	185
	186	hitForRec->SetBendingCoor(y);
	187	hitForRec->SetNonBendingCoor(x);
	188	hitForRec->SetZ(z);
	189	hitForRec->SetBendingReso2(0.0);
	190	hitForRec->SetNonBendingReso2(0.0);
	191	iChamber = AliMUONConstants::ChamberNumber(z);
	192	hitForRec->SetChamberNumber(iChamber);
	193
	194	muonTrack->AddTrackParamAtHit(trackParam);
	195	muonTrack->AddHitForRecAtHit(hitForRec);
	196	muonTrack->SetTrackID(track);
	197
	198	trackSave = track;
	199	if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
	200	}
	201
	202	// track parameters at vertex
	203	particle = fRunLoader->GetHeader()->Stack()->Particle(muonTrack->GetTrackID());
	204	if (particle) {
	205	x = particle->Vx();
	206	y = particle->Vy();
	207	z = particle->Vz();
	208	pX = particle->Px();
	209	pY = particle->Py();
	210	pZ = particle->Pz();
	211
	212	trackParam->SetBendingCoor(y);
	213	trackParam->SetNonBendingCoor(x);
	214	trackParam->SetZ(z);
	215	if (TMath::Abs(pZ) > 0) {
	216	bendingSlope = pY/pZ;
	217	nonBendingSlope = pX/pZ;
	218	}
	219	pYZ = TMath::Sqrt(pYpY+pZpZ);
	220	if (pYZ >0) inverseBendingMomentum = 1/pYZ;
	221	trackParam->SetBendingSlope(bendingSlope);
	222	trackParam->SetNonBendingSlope(nonBendingSlope);
	223	trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
	224
	225	muonTrack->SetTrackParamAtVertex(trackParam);
	226	}
	227
	228	AddMuonTrackReference(muonTrack);
	229	muonTrack->ResetTrackParamAtHit();
	230	muonTrack->ResetHitForRecAtHit();
	231
	232	} // end while isNewTrack
	233
	234	}
	235
	236	CleanMuonTrackRef();
	237
	238	ReconstructibleTracks();
	239
	240	delete muonTrack;
	241	delete trackParam;
	242	delete hitForRec;
	243
	244	}
	245
	246	//____________________________________________________________________________
	247	TClonesArray* AliMUONRecoCheck::GetTrackReco()
	248	{
	249	/// Return TClonesArray of reconstructed tracks
	250
	251	GetMUONData()->ResetRecTracks();
	252	GetMUONData()->SetTreeAddress("RT");
	253	fTrackReco = GetMUONData()->RecTracks();
	254	GetMUONData()->GetRecTracks();
	255	fRecoTracks = fTrackReco->GetEntriesFast();
	256	return fTrackReco;
	257	}
	258	//_____________________________________________________________________________
	259	void AliMUONRecoCheck::PrintEvent() const
	260	{
	261	/// Debug facility
	262
	263	AliMUONTrack *track;
	264	AliMUONHitForRec *hitForRec;
	265	TClonesArray * hitForRecAtHit = 0;
	266	Float_t xRec,yRec,zRec;
	267
	268	Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
	269
	270	printf(" ******************************* \n");
	271	printf(" nb of tracks %d \n",nTrackRef);
	272
	273	for (Int_t index = 0; index < nTrackRef; index++) {
	274	track = (AliMUONTrack*)fMuonTrackRef->At(index);
	275	hitForRecAtHit = track->GetHitForRecAtHit();
	276	Int_t nTrackHits = hitForRecAtHit->GetEntriesFast();
	277	printf(" track number %d \n",index);
	278	for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
	279	hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
	280	xRec = hitForRec->GetNonBendingCoor();
	281	yRec = hitForRec->GetBendingCoor();
	282	zRec = hitForRec->GetZ();
	283	printf(" x,y,z: %f , %f , %f \n",xRec,yRec,zRec);
	284	}
	285	}
	286	}
	287
	288	//_____________________________________________________________________________
	289	void AliMUONRecoCheck::ResetTracks() const
	290	{
	291	/// Reset tracks
	292
	293	if (fMuonTrackRef) fMuonTrackRef->Clear();
	294	}
	295	//_____________________________________________________________________________
	296	void AliMUONRecoCheck::CleanMuonTrackRef()
	297	{
	298	/// Re-calculate hits parameters because two AliTrackReferences are recorded for
	299	/// each chamber (one when particle is entering + one when particle is leaving
	300	/// the sensitive volume)
	301
	302	Float_t maxGasGap = 1.; // cm
	303	AliMUONTrack track, trackNew;
	304	AliMUONHitForRec hitForRec, hitForRec1, *hitForRec2;
	305	AliMUONTrackParam trackParam, trackParam1, trackParam2, trackParamAtVertex;
	306	TClonesArray * hitForRecAtHit = 0;
	307	TClonesArray * trackParamAtHit = 0;
	308	Float_t xRec,yRec,zRec;
	309	Float_t xRec1,yRec1,zRec1;
	310	Float_t xRec2,yRec2,zRec2;
	311	Float_t bendingSlope,nonBendingSlope,bendingMomentum;
	312	Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
	313	Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
	314	TClonesArray *newMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
	315	Int_t iHit1;
	316	Int_t iChamber = 0;
	317	Int_t nRec = 0;
	318	Int_t nTrackHits = 0;
	319
	320	hitForRec = new AliMUONHitForRec();
	321	trackParam = new AliMUONTrackParam();
	322	trackNew = new AliMUONTrack();
	323
	324	Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
	325
	326	for (Int_t index = 0; index < nTrackRef; index++) {
	327	track = (AliMUONTrack*)fMuonTrackRef->At(index);
	328	hitForRecAtHit = track->GetHitForRecAtHit();
	329	trackParamAtHit = track->GetTrackParamAtHit();
	330	trackParamAtVertex = track->GetTrackParamAtVertex();
	331	nTrackHits = hitForRecAtHit->GetEntriesFast();
	332	iHit1 = 0;
	333	while (iHit1 < nTrackHits) {
	334	hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1);
	335	trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1);
	336	xRec1 = hitForRec1->GetNonBendingCoor();
	337	yRec1 = hitForRec1->GetBendingCoor();
	338	zRec1 = hitForRec1->GetZ();
	339	xRec = xRec1;
	340	yRec = yRec1;
	341	zRec = zRec1;
	342	bendingSlope1 = trackParam1->GetBendingSlope();
	343	nonBendingSlope1 = trackParam1->GetNonBendingSlope();
	344	bendingMomentum1 = 0;
	345	if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
	346	bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
	347	bendingSlope = bendingSlope1;
	348	nonBendingSlope = nonBendingSlope1;
	349	bendingMomentum = bendingMomentum1;
	350	nRec = 1;
	351	for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
	352	hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2);
	353	trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2);
	354	xRec2 = hitForRec2->GetNonBendingCoor();
	355	yRec2 = hitForRec2->GetBendingCoor();
	356	zRec2 = hitForRec2->GetZ();
	357	bendingSlope2 = trackParam2->GetBendingSlope();
	358	nonBendingSlope2 = trackParam2->GetNonBendingSlope();
	359	bendingMomentum2 = 0;
	360	if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
	361	bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
	362
	363	if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
	364	nRec++;
	365	xRec += xRec2;
	366	yRec += yRec2;
	367	zRec += zRec2;
	368	bendingSlope += bendingSlope2;
	369	nonBendingSlope += nonBendingSlope2;
	370	bendingMomentum += bendingMomentum2;
	371	iHit1 = iHit2;
	372	}
	373
	374	} // end iHit2
	375	xRec /= (Float_t)nRec;
	376	yRec /= (Float_t)nRec;
	377	zRec /= (Float_t)nRec;
	378	bendingSlope /= (Float_t)nRec;
	379	nonBendingSlope /= (Float_t)nRec;
	380	bendingMomentum /= (Float_t)nRec;
	381
	382	hitForRec->SetNonBendingCoor(xRec);
	383	hitForRec->SetBendingCoor(yRec);
	384	hitForRec->SetZ(zRec);
	385	iChamber = AliMUONConstants::ChamberNumber(zRec);
	386	hitForRec->SetChamberNumber(iChamber);
	387	hitForRec->SetBendingReso2(0.0);
	388	hitForRec->SetNonBendingReso2(0.0);
	389	trackParam->SetNonBendingCoor(xRec);
	390	trackParam->SetBendingCoor(yRec);
	391	trackParam->SetZ(zRec);
	392	trackParam->SetNonBendingSlope(nonBendingSlope);
	393	trackParam->SetBendingSlope(bendingSlope);
	394	if (TMath::Abs(bendingMomentum) > 0)
	395	trackParam->SetInverseBendingMomentum(1./bendingMomentum);
	396
	397	trackNew->AddHitForRecAtHit(hitForRec);
	398	trackNew->AddTrackParamAtHit(trackParam);
	399
	400	iHit1++;
	401	} // end iHit1
	402
	403	trackNew->SetTrackID(track->GetTrackID());
	404	trackNew->SetTrackParamAtVertex(trackParamAtVertex);
	405	{new ((newMuonTrackRef)[newMuonTrackRef->GetEntriesFast()]) AliMUONTrack(trackNew);}
	406	trackNew->ResetHitForRecAtHit();
	407	trackNew->ResetTrackParamAtHit();
	408
	409	} // end trackRef
	410
	411	fMuonTrackRef->Clear();
	412	nTrackRef = newMuonTrackRef->GetEntriesFast();
	413	for (Int_t index = 0; index < nTrackRef; index++) {
	414	track = (AliMUONTrack*)newMuonTrackRef->At(index);
	415	AddMuonTrackReference(track);
	416	}
	417
	418
	419	delete trackNew;
	420	delete hitForRec;
	421	delete trackParam;
	422	newMuonTrackRef->Delete();
	423	delete newMuonTrackRef;
	424
	425	}
	426
	427	//_____________________________________________________________________________
	428	void AliMUONRecoCheck::ReconstructibleTracks()
	429	{
	430	/// Calculate the number of reconstructible tracks
	431
	432	AliMUONTrack* track;
	433	TClonesArray* hitForRecAtHit = NULL;
	434	AliMUONHitForRec* hitForRec;
	435	Float_t zRec;
	436	Int_t nTrackRef, nTrackHits;
	437	Int_t isChamberInTrack[10];
	438	Int_t iChamber = 0;
	439	Bool_t isTrackOK = kTRUE;
	440
	441	fReconstructibleTracks = 0;
	442
	443	nTrackRef = fMuonTrackRef->GetEntriesFast();
	444	for (Int_t index = 0; index < nTrackRef; index++) {
	445	track = (AliMUONTrack*)fMuonTrackRef->At(index);
	446	hitForRecAtHit = track->GetHitForRecAtHit();
	447	nTrackHits = hitForRecAtHit->GetEntriesFast();
	448	for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
	449
	450	for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
	451	hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
	452	zRec = hitForRec->GetZ();
	453	iChamber = hitForRec->GetChamberNumber();
	454	if (iChamber < 0 \|\| iChamber > 10) continue;
	455	isChamberInTrack[iChamber] = 1;
	456
	457	}
	458	// track is reconstructible if the particle is crossing every tracking chambers
	459	isTrackOK = kTRUE;
	460	for (Int_t ch = 0; ch < 10; ch++) {
	461	if (!isChamberInTrack[ch]) isTrackOK = kFALSE;
	462	} if (isTrackOK) fReconstructibleTracks++;
	463	if (!isTrackOK) fMuonTrackRef->Remove(track); // remove non reconstructible tracks
	464	}
	465	fMuonTrackRef->Compress();
	466	}
	467
	468