[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
/// \brief Utility class to check 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 "AliMUON.h"
#include "AliMUONRecoCheck.h"
#include "AliMUONTrack.h"
#include "AliMUONData.h"
#include "AliMUONConstants.h"

#include "AliRun.h" // for gAlice
#include "AliLoader.h" 
#include "AliRunLoader.h" 
#include "AliHeader.h"
#include "AliStack.h"
#include "AliTrackReference.h"
#include "AliLog.h" 

#include <TParticle.h>

/// \cond CLASSIMP
ClassImp(AliMUONRecoCheck)
/// \endcond

//_____________________________________________________________________________
  AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader)
  : TObject(),
  fRunLoader(0x0),
  fMUONData(0x0),
  fMuonTrackRef(0x0),
  fTrackReco(0x0),
  fReconstructibleTracks(0),
  fRecoTracks(0)
{
/// 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");

}

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

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

//_____________________________________________________________________________
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, detElemId;

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

  Int_t max = fRunLoader->GetHeader()->Stack()->GetNtrack();
  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 >= max){
	  AliWarningStream()
	    << "Track ID " << track 
	    << " larger than max number of particles " << max << endl;
	  isNewTrack = kFALSE;
	  break;
	}
	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);
	detElemId = hitForRec->GetDetElemId();
	if (detElemId) iChamber = detElemId / 100 - 1; 
	else 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, detElemId = 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);
      detElemId = hitForRec->GetDetElemId();
      if (detElemId) iChamber = detElemId / 100 - 1;
      else 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 depositing a hit
    // in the following chamber combinations:

    isTrackOK = kTRUE;
    if (!isChamberInTrack[0] && !isChamberInTrack[1]) isTrackOK = kFALSE;
    if (!isChamberInTrack[2] && !isChamberInTrack[3]) isTrackOK = kFALSE;
    if (!isChamberInTrack[4] && !isChamberInTrack[5]) isTrackOK = kFALSE;
    Int_t nHitsInLastStations=0;
    for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++)
      if (isChamberInTrack[ch]) nHitsInLastStations++; 
    if(nHitsInLastStations < 3) isTrackOK = kFALSE;

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