[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 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 "AliMUONRecoCheck.h"
#include "AliMUONHitForRec.h"
#include "AliMUONTrack.h"
#include "AliMUONConstants.h"
#include "AliMUONMCDataInterface.h"
#include "AliMUONDataInterface.h"
#include "AliStack.h"
#include "AliTrackReference.h"
#include "AliLog.h" 
#include "AliMUONTrackStoreV1.h"
#include <Riostream.h>
#include <TParticle.h>
#include <TParticlePDG.h>

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

//_____________________________________________________________________________
AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader, Char_t *chLoaderSim)
: TObject(),
fMCDataInterface(new AliMUONMCDataInterface(chLoaderSim)),
fDataInterface(new AliMUONDataInterface(chLoader))
{
  /// Normal ctor
}

//_____________________________________________________________________________
AliMUONRecoCheck::~AliMUONRecoCheck()
{
  /// Destructor
  delete fMCDataInterface;
  delete fDataInterface;
}

//_____________________________________________________________________________
AliMUONVTrackStore* 
AliMUONRecoCheck::ReconstructedTracks(Int_t event)
{
  /// Return the reconstructed track store for a given event
  return fDataInterface->TrackStore(event);
}

//_____________________________________________________________________________
AliMUONVTrackStore* 
AliMUONRecoCheck::TrackRefs(Int_t event)
{
  /// Return a track store containing the track references (converted into 
  /// MUONTrack objects) for a given event
  return MakeTrackRefs(event);
}

//_____________________________________________________________________________
AliMUONVTrackStore* 
AliMUONRecoCheck::ReconstructibleTracks(Int_t event)
{
  /// Return a track store containing the reconstructible tracks for a given event
  AliMUONVTrackStore* tmp = MakeTrackRefs(event);
  AliMUONVTrackStore* reconstructible = MakeReconstructibleTracks(*tmp);
  delete tmp;
  return reconstructible;
}

//_____________________________________________________________________________
AliMUONVTrackStore*
AliMUONRecoCheck::MakeTrackRefs(Int_t event)
{
  /// Make reconstructible tracks
    
  AliMUONVTrackStore* trackRefStore = new AliMUONTrackStoreV1;
  
  Int_t nTrackRef = fMCDataInterface->NumberOfTrackRefs(event);
  
  Int_t trackSave(-999);
  Int_t iHitMin(0);
  
  Bool_t isNewTrack;
    
  AliStack* stack = fMCDataInterface->Stack(event);
  Int_t max = stack->GetNtrack();
  
  for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; ++iTrackRef) 
  {
    TClonesArray* trackRefs = fMCDataInterface->TrackRefs(event,iTrackRef);
    
    iHitMin = 0;
    isNewTrack = kTRUE;
    
    if (!trackRefs->GetEntries()) continue; 
    
    while (isNewTrack) {
      
      AliMUONTrack muonTrack;
      
      for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); ++iHit) 
      {        
        AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->At(iHit));
        
        Float_t x = trackReference->X();
        Float_t y = trackReference->Y();
        Float_t z = trackReference->Z();
        Float_t pX = trackReference->Px();
        Float_t pY = trackReference->Py();
        Float_t pZ = trackReference->Pz();
        
        Int_t 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;
        }
        
        Float_t bendingSlope = 0;
        Float_t nonBendingSlope = 0;
        Float_t inverseBendingMomentum = 0;
        
        AliMUONTrackParam trackParam;
        
        // track parameters at hit
        trackParam.SetBendingCoor(y);
        trackParam.SetNonBendingCoor(x);
        trackParam.SetZ(z);
        
        if (TMath::Abs(pZ) > 0) 
        {
          bendingSlope = pY/pZ;
          nonBendingSlope = pX/pZ;
        }
        Float_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
        if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
        
        trackParam.SetBendingSlope(bendingSlope);
        trackParam.SetNonBendingSlope(nonBendingSlope);
        trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
        
        AliMUONHitForRec hitForRec;
        
        hitForRec.SetBendingCoor(y);
        hitForRec.SetNonBendingCoor(x);
        hitForRec.SetZ(z);
        hitForRec.SetBendingReso2(0.0); 
        hitForRec.SetNonBendingReso2(0.0);
        Int_t detElemId = hitForRec.GetDetElemId();
        Int_t iChamber;
        if (detElemId) 
        {
          iChamber = detElemId / 100 - 1; 
        }
        else 
        {
          iChamber = AliMUONConstants::ChamberNumber(z);
        }
        hitForRec.SetChamberNumber(iChamber);
        
        muonTrack.AddTrackParamAtHit(&trackParam,0);
        muonTrack.AddHitForRecAtHit(&hitForRec);
        muonTrack.SetTrackID(track);
        
        trackSave = track;
        if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
      }
      
      // track parameters at vertex 
      TParticle* particle = stack->Particle(muonTrack.GetTrackID());
      
      if (particle) 
      {        
        Float_t x = particle->Vx();
        Float_t y = particle->Vy();
        Float_t z = particle->Vz();
        Float_t pX = particle->Px();
        Float_t pY = particle->Py();
        Float_t pZ = particle->Pz();
        
        AliMUONTrackParam trackParam;
        
        trackParam.SetBendingCoor(y);
        trackParam.SetNonBendingCoor(x);
        trackParam.SetZ(z);
        
        Float_t bendingSlope = 0;
        Float_t nonBendingSlope = 0;
        Float_t inverseBendingMomentum = 0;
                
        if (TMath::Abs(pZ) > 0) 
        {
          bendingSlope = pY/pZ;
          nonBendingSlope = pX/pZ;
        }
        
        Float_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
        if (pYZ >0) inverseBendingMomentum = 1/pYZ;       
        
        TParticlePDG* ppdg = particle->GetPDG(1);
        Int_t charge = (Int_t)(ppdg->Charge()/3.0);
        inverseBendingMomentum *= charge;
        
        trackParam.SetBendingSlope(bendingSlope);
        trackParam.SetNonBendingSlope(nonBendingSlope);
        trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
        
        muonTrack.SetTrackParamAtVertex(&trackParam);
      }

      trackRefStore->Add(muonTrack);
    } // end while isNewTrack
    delete trackRefs;
  }
  
  AliMUONVTrackStore* rv = CleanMuonTrackRef(*trackRefStore);
  
  delete trackRefStore;

  return rv;
}

//_____________________________________________________________________________
Int_t
AliMUONRecoCheck::NumberOfEvents() const
{
  /// Return the number of events
  if ( fDataInterface ) 
  {
    return fDataInterface->NumberOfEvents();
  }
  return 0;
}

//_____________________________________________________________________________
AliMUONVTrackStore*
AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore& trackRefs)
{
  /// 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 
  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;
  
  AliMUONVTrackStore* newMuonTrackRef = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
  Int_t iHit1;
  Int_t iChamber = 0, detElemId = 0;
  Int_t nRec = 0;
  Int_t nTrackHits = 0;
  
  hitForRec = new AliMUONHitForRec();
  trackParam = new AliMUONTrackParam();
  
  TIter next(trackRefs.CreateIterator());
  AliMUONTrack* track;
  
  while ( ( track = static_cast<AliMUONTrack*>(next())) ) 
  {
    hitForRecAtHit = track->GetHitForRecAtHit();
    trackParamAtHit = track->GetTrackParamAtHit();
    trackParamAtVertex = track->GetTrackParamAtVertex();
    nTrackHits = hitForRecAtHit->GetEntriesFast();
    AliMUONTrack trackNew;
    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,0);
      
      iHit1++;
    } // end iHit1
    
    trackNew.SetTrackID(track->GetTrackID());
    trackNew.SetTrackParamAtVertex(trackParamAtVertex);
    newMuonTrackRef->Add(trackNew);
    
  } // end trackRef
  
  delete hitForRec;
  delete trackParam;
  return newMuonTrackRef;  
}

//_____________________________________________________________________________
AliMUONVTrackStore*
AliMUONRecoCheck::MakeReconstructibleTracks(const AliMUONVTrackStore& trackRefs)
{
  /// Calculate the number of reconstructible tracks
  
  AliMUONVTrackStore* reconstructibleStore = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
  
  TClonesArray* hitForRecAtHit = NULL;
  AliMUONHitForRec* hitForRec;
  Float_t zRec;
  Int_t nTrackHits;
  Int_t isChamberInTrack[10];
  Int_t iChamber = 0;
  Bool_t isTrackOK = kTRUE;
    
  TIter next(trackRefs.CreateIterator());
  AliMUONTrack* track;

  while ( ( track = static_cast<AliMUONTrack*>(next()) ) )
  {
    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) 
    {
      reconstructibleStore->Add(*track);
    }
  }

  return reconstructibleStore;
}
Commit	Line	Data
b8dc484b	1	/**************************************************************************
48217459	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	**************************************************************************/
b8dc484b	15
d19b6003	16	/* $Id$ */
d19b6003	17
3d1463c8	18	//-----------------------------------------------------------------------------
5398f946	19	/// \class AliMUONRecoCheck
78649106	20	/// Utility class to check reconstruction
5398f946	21	/// Reconstructed tracks are compared to reference tracks.
	22	/// The reference tracks are built from AliTrackReference for the
	23	/// hit in chamber (0..9) and from kinematics for the vertex parameters.
3d1463c8	24	//-----------------------------------------------------------------------------
b8dc484b	25
cc056de1	26	#include "AliMUONRecoCheck.h"
b16234ab	27	#include "AliMUONHitForRec.h"
cc056de1	28	#include "AliMUONTrack.h"
cc056de1	29	#include "AliMUONConstants.h"
48217459	30	#include "AliMUONMCDataInterface.h"
48217459	31	#include "AliMUONDataInterface.h"
b8dc484b	32	#include "AliStack.h"
cc056de1	33	#include "AliTrackReference.h"
cc056de1	34	#include "AliLog.h"
48217459	35	#include "AliMUONTrackStoreV1.h"
48217459	36	#include <Riostream.h>
cc056de1	37	#include <TParticle.h>
5b1afea0	38	#include <TParticlePDG.h>
b8dc484b	39
5398f946	40	/// \cond CLASSIMP
b8dc484b	41	ClassImp(AliMUONRecoCheck)
5398f946	42	/// \endcond
b8dc484b	43
068e4c36	44	//_____________________________________________________________________________
48217459	45	AliMUONRecoCheck::AliMUONRecoCheck(Char_t chLoader, Char_t chLoaderSim)
	46	: TObject(),
	47	fMCDataInterface(new AliMUONMCDataInterface(chLoaderSim)),
	48	fDataInterface(new AliMUONDataInterface(chLoader))
068e4c36	49	{
48217459	50	/// Normal ctor
068e4c36	51	}
068e4c36	52
b8dc484b	53	//_____________________________________________________________________________
48217459	54	AliMUONRecoCheck::~AliMUONRecoCheck()
b8dc484b	55	{
48217459	56	/// Destructor
	57	delete fMCDataInterface;
	58	delete fDataInterface;
b8dc484b	59	}
b8dc484b	60
48217459	61	//_____________________________________________________________________________
	62	AliMUONVTrackStore*
	63	AliMUONRecoCheck::ReconstructedTracks(Int_t event)
	64	{
	65	/// Return the reconstructed track store for a given event
	66	return fDataInterface->TrackStore(event);
	67	}
45a05f07	68
b8dc484b	69	//_____________________________________________________________________________
48217459	70	AliMUONVTrackStore*
48217459	71	AliMUONRecoCheck::TrackRefs(Int_t event)
b8dc484b	72	{
48217459	73	/// Return a track store containing the track references (converted into
	74	/// MUONTrack objects) for a given event
	75	return MakeTrackRefs(event);
b8dc484b	76	}
	77
	78	//_____________________________________________________________________________
48217459	79	AliMUONVTrackStore*
48217459	80	AliMUONRecoCheck::ReconstructibleTracks(Int_t event)
b8dc484b	81	{
48217459	82	/// Return a track store containing the reconstructible tracks for a given event
	83	AliMUONVTrackStore* tmp = MakeTrackRefs(event);
	84	AliMUONVTrackStore* reconstructible = MakeReconstructibleTracks(*tmp);
	85	delete tmp;
	86	return reconstructible;
	87	}
b8dc484b	88
48217459	89	//_____________________________________________________________________________
	90	AliMUONVTrackStore*
	91	AliMUONRecoCheck::MakeTrackRefs(Int_t event)
	92	{
	93	/// Make reconstructible tracks
	94
	95	AliMUONVTrackStore* trackRefStore = new AliMUONTrackStoreV1;
	96
	97	Int_t nTrackRef = fMCDataInterface->NumberOfTrackRefs(event);
	98
	99	Int_t trackSave(-999);
	100	Int_t iHitMin(0);
b8dc484b	101
b8dc484b	102	Bool_t isNewTrack;
48217459	103
	104	AliStack* stack = fMCDataInterface->Stack(event);
	105	Int_t max = stack->GetNtrack();
	106
	107	for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; ++iTrackRef)
	108	{
	109	TClonesArray* trackRefs = fMCDataInterface->TrackRefs(event,iTrackRef);
b8dc484b	110
	111	iHitMin = 0;
	112	isNewTrack = kTRUE;
dba3c809	113
29fc2c86	114	if (!trackRefs->GetEntries()) continue;
48217459	115
29fc2c86	116	while (isNewTrack) {
068e4c36	117
48217459	118	AliMUONTrack muonTrack;
b8dc484b	119
48217459	120	for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); ++iHit)
	121	{
	122	AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->At(iHit));
	123
	124	Float_t x = trackReference->X();
	125	Float_t y = trackReference->Y();
	126	Float_t z = trackReference->Z();
	127	Float_t pX = trackReference->Px();
	128	Float_t pY = trackReference->Py();
	129	Float_t pZ = trackReference->Pz();
	130
	131	Int_t track = trackReference->GetTrack();
	132
	133	if ( track >= max )
	134	{
	135	AliWarningStream()
	136	<< "Track ID " << track
	137	<< " larger than max number of particles " << max << endl;
	138	isNewTrack = kFALSE;
	139	break;
	140	}
	141	if (track != trackSave && iHit != 0) {
	142	iHitMin = iHit;
	143	trackSave = track;
	144	break;
	145	}
	146
	147	Float_t bendingSlope = 0;
	148	Float_t nonBendingSlope = 0;
	149	Float_t inverseBendingMomentum = 0;
	150
	151	AliMUONTrackParam trackParam;
	152
	153	// track parameters at hit
	154	trackParam.SetBendingCoor(y);
	155	trackParam.SetNonBendingCoor(x);
	156	trackParam.SetZ(z);
	157
	158	if (TMath::Abs(pZ) > 0)
	159	{
	160	bendingSlope = pY/pZ;
	161	nonBendingSlope = pX/pZ;
	162	}
	163	Float_t pYZ = TMath::Sqrt(pYpY+pZpZ);
	164	if (pYZ >0) inverseBendingMomentum = 1/pYZ;
	165
	166	trackParam.SetBendingSlope(bendingSlope);
	167	trackParam.SetNonBendingSlope(nonBendingSlope);
	168	trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
	169
	170	AliMUONHitForRec hitForRec;
	171
	172	hitForRec.SetBendingCoor(y);
	173	hitForRec.SetNonBendingCoor(x);
	174	hitForRec.SetZ(z);
	175	hitForRec.SetBendingReso2(0.0);
	176	hitForRec.SetNonBendingReso2(0.0);
	177	Int_t detElemId = hitForRec.GetDetElemId();
	178	Int_t iChamber;
	179	if (detElemId)
	180	{
	181	iChamber = detElemId / 100 - 1;
	182	}
	183	else
184	{
185	iChamber = AliMUONConstants::ChamberNumber(z);
186	}
187	hitForRec.SetChamberNumber(iChamber);
188
189	muonTrack.AddTrackParamAtHit(&trackParam,0);
190	muonTrack.AddHitForRecAtHit(&hitForRec);
191	muonTrack.SetTrackID(track);
192
193	trackSave = track;
194	if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
b8dc484b	195	}
48217459	196
b8dc484b	197	// track parameters at vertex
48217459	198	TParticle* particle = stack->Particle(muonTrack.GetTrackID());
b8dc484b	199
48217459	200	if (particle)
	201	{
	202	Float_t x = particle->Vx();
	203	Float_t y = particle->Vy();
	204	Float_t z = particle->Vz();
	205	Float_t pX = particle->Px();
	206	Float_t pY = particle->Py();
	207	Float_t pZ = particle->Pz();
	208
	209	AliMUONTrackParam trackParam;
	210
	211	trackParam.SetBendingCoor(y);
	212	trackParam.SetNonBendingCoor(x);
	213	trackParam.SetZ(z);
	214
	215	Float_t bendingSlope = 0;
	216	Float_t nonBendingSlope = 0;
	217	Float_t inverseBendingMomentum = 0;
	218
	219	if (TMath::Abs(pZ) > 0)
	220	{
	221	bendingSlope = pY/pZ;
	222	nonBendingSlope = pX/pZ;
	223	}
	224
	225	Float_t pYZ = TMath::Sqrt(pYpY+pZpZ);
	226	if (pYZ >0) inverseBendingMomentum = 1/pYZ;
	227
	228	TParticlePDG* ppdg = particle->GetPDG(1);
	229	Int_t charge = (Int_t)(ppdg->Charge()/3.0);
	230	inverseBendingMomentum *= charge;
	231
	232	trackParam.SetBendingSlope(bendingSlope);
	233	trackParam.SetNonBendingSlope(nonBendingSlope);
	234	trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
	235
	236	muonTrack.SetTrackParamAtVertex(&trackParam);
b8dc484b	237	}
cc056de1	238
48217459	239	trackRefStore->Add(muonTrack);
b8dc484b	240	} // end while isNewTrack
48217459	241	delete trackRefs;
b8dc484b	242	}
b8dc484b	243
48217459	244	AliMUONVTrackStore* rv = CleanMuonTrackRef(*trackRefStore);
	245
	246	delete trackRefStore;
b8dc484b	247
48217459	248	return rv;
b8dc484b	249	}
b8dc484b	250
b8dc484b	251	//_____________________________________________________________________________
48217459	252	Int_t
48217459	253	AliMUONRecoCheck::NumberOfEvents() const
b8dc484b	254	{
48217459	255	/// Return the number of events
	256	if ( fDataInterface )
	257	{
	258	return fDataInterface->NumberOfEvents();
b8dc484b	259	}
48217459	260	return 0;
b8dc484b	261	}
	262
	263	//_____________________________________________________________________________
48217459	264	AliMUONVTrackStore*
48217459	265	AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore& trackRefs)
b8dc484b	266	{
48217459	267	/// Re-calculate hits parameters because two AliTrackReferences are recorded for
	268	/// each chamber (one when particle is entering + one when particle is leaving
	269	/// the sensitive volume)
	270
b8dc484b	271	Float_t maxGasGap = 1.; // cm
b8dc484b	272	AliMUONHitForRec hitForRec, hitForRec1, *hitForRec2;
	273	AliMUONTrackParam trackParam, trackParam1, trackParam2, trackParamAtVertex;
	274	TClonesArray * hitForRecAtHit = 0;
	275	TClonesArray * trackParamAtHit = 0;
	276	Float_t xRec,yRec,zRec;
	277	Float_t xRec1,yRec1,zRec1;
	278	Float_t xRec2,yRec2,zRec2;
	279	Float_t bendingSlope,nonBendingSlope,bendingMomentum;
	280	Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
	281	Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
48217459	282
48217459	283	AliMUONVTrackStore* newMuonTrackRef = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
b8dc484b	284	Int_t iHit1;
30d3ea8a	285	Int_t iChamber = 0, detElemId = 0;
b8dc484b	286	Int_t nRec = 0;
b8dc484b	287	Int_t nTrackHits = 0;
48217459	288
b8dc484b	289	hitForRec = new AliMUONHitForRec();
b8dc484b	290	trackParam = new AliMUONTrackParam();
48217459	291
	292	TIter next(trackRefs.CreateIterator());
	293	AliMUONTrack* track;
	294
	295	while ( ( track = static_cast<AliMUONTrack*>(next())) )
	296	{
b8dc484b	297	hitForRecAtHit = track->GetHitForRecAtHit();
	298	trackParamAtHit = track->GetTrackParamAtHit();
	299	trackParamAtVertex = track->GetTrackParamAtVertex();
	300	nTrackHits = hitForRecAtHit->GetEntriesFast();
48217459	301	AliMUONTrack trackNew;
b8dc484b	302	iHit1 = 0;
48217459	303	while (iHit1 < nTrackHits)
48217459	304	{
b8dc484b	305	hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1);
	306	trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1);
	307	xRec1 = hitForRec1->GetNonBendingCoor();
	308	yRec1 = hitForRec1->GetBendingCoor();
	309	zRec1 = hitForRec1->GetZ();
	310	xRec = xRec1;
	311	yRec = yRec1;
	312	zRec = zRec1;
	313	bendingSlope1 = trackParam1->GetBendingSlope();
	314	nonBendingSlope1 = trackParam1->GetNonBendingSlope();
de876dd1	315	bendingMomentum1 = 0;
de876dd1	316	if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
48217459	317	bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
b8dc484b	318	bendingSlope = bendingSlope1;
	319	nonBendingSlope = nonBendingSlope1;
	320	bendingMomentum = bendingMomentum1;
	321	nRec = 1;
48217459	322	for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++)
	323	{
	324	hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2);
	325	trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2);
	326	xRec2 = hitForRec2->GetNonBendingCoor();
	327	yRec2 = hitForRec2->GetBendingCoor();
	328	zRec2 = hitForRec2->GetZ();
	329	bendingSlope2 = trackParam2->GetBendingSlope();
	330	nonBendingSlope2 = trackParam2->GetNonBendingSlope();
	331	bendingMomentum2 = 0;
	332	if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
	333	bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
	334
	335	if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
	336
	337	nRec++;
	338	xRec += xRec2;
	339	yRec += yRec2;
	340	zRec += zRec2;
	341	bendingSlope += bendingSlope2;
	342	nonBendingSlope += nonBendingSlope2;
	343	bendingMomentum += bendingMomentum2;
	344	iHit1 = iHit2;
	345	}
	346
b8dc484b	347	} // end iHit2
	348	xRec /= (Float_t)nRec;
	349	yRec /= (Float_t)nRec;
	350	zRec /= (Float_t)nRec;
	351	bendingSlope /= (Float_t)nRec;
	352	nonBendingSlope /= (Float_t)nRec;
	353	bendingMomentum /= (Float_t)nRec;
	354
	355	hitForRec->SetNonBendingCoor(xRec);
	356	hitForRec->SetBendingCoor(yRec);
	357	hitForRec->SetZ(zRec);
30d3ea8a	358	detElemId = hitForRec->GetDetElemId();
	359	if (detElemId) iChamber = detElemId / 100 - 1;
	360	else iChamber = AliMUONConstants::ChamberNumber(zRec);
b8dc484b	361	hitForRec->SetChamberNumber(iChamber);
	362	hitForRec->SetBendingReso2(0.0);
	363	hitForRec->SetNonBendingReso2(0.0);
	364	trackParam->SetNonBendingCoor(xRec);
	365	trackParam->SetBendingCoor(yRec);
	366	trackParam->SetZ(zRec);
	367	trackParam->SetNonBendingSlope(nonBendingSlope);
	368	trackParam->SetBendingSlope(bendingSlope);
de876dd1	369	if (TMath::Abs(bendingMomentum) > 0)
48217459	370	trackParam->SetInverseBendingMomentum(1./bendingMomentum);
	371
	372	trackNew.AddHitForRecAtHit(hitForRec);
	373	trackNew.AddTrackParamAtHit(trackParam,0);
b8dc484b	374
	375	iHit1++;
	376	} // end iHit1
48217459	377
	378	trackNew.SetTrackID(track->GetTrackID());
	379	trackNew.SetTrackParamAtVertex(trackParamAtVertex);
	380	newMuonTrackRef->Add(trackNew);
b8dc484b	381
b8dc484b	382	} // end trackRef
b8dc484b	383
b8dc484b	384	delete hitForRec;
b8dc484b	385	delete trackParam;
48217459	386	return newMuonTrackRef;
b8dc484b	387	}
	388
	389	//_____________________________________________________________________________
48217459	390	AliMUONVTrackStore*
48217459	391	AliMUONRecoCheck::MakeReconstructibleTracks(const AliMUONVTrackStore& trackRefs)
b8dc484b	392	{
48217459	393	/// Calculate the number of reconstructible tracks
	394
	395	AliMUONVTrackStore* reconstructibleStore = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
b8dc484b	396
b8dc484b	397	TClonesArray* hitForRecAtHit = NULL;
	398	AliMUONHitForRec* hitForRec;
	399	Float_t zRec;
48217459	400	Int_t nTrackHits;
b8dc484b	401	Int_t isChamberInTrack[10];
	402	Int_t iChamber = 0;
	403	Bool_t isTrackOK = kTRUE;
48217459	404
	405	TIter next(trackRefs.CreateIterator());
	406	AliMUONTrack* track;
b8dc484b	407
48217459	408	while ( ( track = static_cast<AliMUONTrack*>(next()) ) )
48217459	409	{
b8dc484b	410	hitForRecAtHit = track->GetHitForRecAtHit();
	411	nTrackHits = hitForRecAtHit->GetEntriesFast();
	412	for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
48217459	413
b8dc484b	414	for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
	415	hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
	416	zRec = hitForRec->GetZ();
48217459	417	iChamber = hitForRec->GetChamberNumber();
b8dc484b	418	if (iChamber < 0 \|\| iChamber > 10) continue;
b8dc484b	419	isChamberInTrack[iChamber] = 1;
cc056de1	420	}
	421	// track is reconstructible if the particle is depositing a hit
	422	// in the following chamber combinations:
48217459	423
b8dc484b	424	isTrackOK = kTRUE;
cc056de1	425	if (!isChamberInTrack[0] && !isChamberInTrack[1]) isTrackOK = kFALSE;
	426	if (!isChamberInTrack[2] && !isChamberInTrack[3]) isTrackOK = kFALSE;
	427	if (!isChamberInTrack[4] && !isChamberInTrack[5]) isTrackOK = kFALSE;
	428	Int_t nHitsInLastStations=0;
	429	for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++)
	430	if (isChamberInTrack[ch]) nHitsInLastStations++;
	431	if(nHitsInLastStations < 3) isTrackOK = kFALSE;
48217459	432
	433	if (isTrackOK)
	434	{
	435	reconstructibleStore->Add(*track);
	436	}
b8dc484b	437	}
48217459	438
48217459	439	return reconstructibleStore;
b8dc484b	440	}
b8dc484b	441