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

//////////////////////////////////////////////////////////////////////////////////////
//                                                                     
// 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 "AliMUONChamber.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()
{
  fRunLoader->UnloadKinematics();
  fRunLoader->UnloadTrackRefs();
  fRunLoader->UnloadTracks();
  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 = new TClonesArray("AliTrackReference", 10);
  branch->SetAddress(&trackRefs);

  Int_t nTrackRef = (Int_t)treeTR->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++) {
    treeTR->GetEntry(iTrackRef);
    
    iHitMin = 0;
    isNewTrack = kTRUE;

    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 = 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;
  delete trackRefs;

}

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

  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
{
  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 = 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 = 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 = 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);
      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;
  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();
}


//_____________________________________________________________________________
Int_t AliMUONRecoCheck::ChamberNumber(Float_t z) const
{
  // return chamber number according z position of hit. Should be taken from geometry ?
 
  Float_t dMaxChamber =  AliMUONConstants::DzSlat() + AliMUONConstants::DzCh() + 0.25; // cm st 3 &4 & 5
  if ( z >  (AliMUONConstants::DefaultChamberZ(4)+50.)) dMaxChamber = 7.; // cm stations 1 & 2
  Int_t iChamber;

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