[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
  }
  
  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
b8dc484b	241	}
b8dc484b	242
48217459	243	AliMUONVTrackStore* rv = CleanMuonTrackRef(*trackRefStore);
	244
	245	delete trackRefStore;
b8dc484b	246
48217459	247	return rv;
b8dc484b	248	}
b8dc484b	249
b8dc484b	250	//_____________________________________________________________________________
48217459	251	Int_t
48217459	252	AliMUONRecoCheck::NumberOfEvents() const
b8dc484b	253	{
48217459	254	/// Return the number of events
	255	if ( fDataInterface )
	256	{
	257	return fDataInterface->NumberOfEvents();
b8dc484b	258	}
48217459	259	return 0;
b8dc484b	260	}
	261
	262	//_____________________________________________________________________________
48217459	263	AliMUONVTrackStore*
48217459	264	AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore& trackRefs)
b8dc484b	265	{
48217459	266	/// Re-calculate hits parameters because two AliTrackReferences are recorded for
	267	/// each chamber (one when particle is entering + one when particle is leaving
	268	/// the sensitive volume)
	269
b8dc484b	270	Float_t maxGasGap = 1.; // cm
b8dc484b	271	AliMUONHitForRec hitForRec, hitForRec1, *hitForRec2;
	272	AliMUONTrackParam trackParam, trackParam1, trackParam2, trackParamAtVertex;
	273	TClonesArray * hitForRecAtHit = 0;
	274	TClonesArray * trackParamAtHit = 0;
	275	Float_t xRec,yRec,zRec;
	276	Float_t xRec1,yRec1,zRec1;
	277	Float_t xRec2,yRec2,zRec2;
	278	Float_t bendingSlope,nonBendingSlope,bendingMomentum;
	279	Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
	280	Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
48217459	281
48217459	282	AliMUONVTrackStore* newMuonTrackRef = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
b8dc484b	283	Int_t iHit1;
30d3ea8a	284	Int_t iChamber = 0, detElemId = 0;
b8dc484b	285	Int_t nRec = 0;
b8dc484b	286	Int_t nTrackHits = 0;
48217459	287
b8dc484b	288	hitForRec = new AliMUONHitForRec();
b8dc484b	289	trackParam = new AliMUONTrackParam();
48217459	290
	291	TIter next(trackRefs.CreateIterator());
	292	AliMUONTrack* track;
	293
	294	while ( ( track = static_cast<AliMUONTrack*>(next())) )
	295	{
b8dc484b	296	hitForRecAtHit = track->GetHitForRecAtHit();
	297	trackParamAtHit = track->GetTrackParamAtHit();
	298	trackParamAtVertex = track->GetTrackParamAtVertex();
	299	nTrackHits = hitForRecAtHit->GetEntriesFast();
48217459	300	AliMUONTrack trackNew;
b8dc484b	301	iHit1 = 0;
48217459	302	while (iHit1 < nTrackHits)
48217459	303	{
b8dc484b	304	hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1);
	305	trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1);
	306	xRec1 = hitForRec1->GetNonBendingCoor();
	307	yRec1 = hitForRec1->GetBendingCoor();
	308	zRec1 = hitForRec1->GetZ();
	309	xRec = xRec1;
	310	yRec = yRec1;
	311	zRec = zRec1;
	312	bendingSlope1 = trackParam1->GetBendingSlope();
	313	nonBendingSlope1 = trackParam1->GetNonBendingSlope();
de876dd1	314	bendingMomentum1 = 0;
de876dd1	315	if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
48217459	316	bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
b8dc484b	317	bendingSlope = bendingSlope1;
	318	nonBendingSlope = nonBendingSlope1;
	319	bendingMomentum = bendingMomentum1;
	320	nRec = 1;
48217459	321	for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++)
	322	{
	323	hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2);
	324	trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2);
	325	xRec2 = hitForRec2->GetNonBendingCoor();
	326	yRec2 = hitForRec2->GetBendingCoor();
	327	zRec2 = hitForRec2->GetZ();
	328	bendingSlope2 = trackParam2->GetBendingSlope();
	329	nonBendingSlope2 = trackParam2->GetNonBendingSlope();
	330	bendingMomentum2 = 0;
	331	if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
	332	bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
	333
	334	if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
	335
	336	nRec++;
	337	xRec += xRec2;
	338	yRec += yRec2;
	339	zRec += zRec2;
	340	bendingSlope += bendingSlope2;
	341	nonBendingSlope += nonBendingSlope2;
	342	bendingMomentum += bendingMomentum2;
	343	iHit1 = iHit2;
	344	}
	345
b8dc484b	346	} // end iHit2
	347	xRec /= (Float_t)nRec;
	348	yRec /= (Float_t)nRec;
	349	zRec /= (Float_t)nRec;
	350	bendingSlope /= (Float_t)nRec;
	351	nonBendingSlope /= (Float_t)nRec;
	352	bendingMomentum /= (Float_t)nRec;
	353
	354	hitForRec->SetNonBendingCoor(xRec);
	355	hitForRec->SetBendingCoor(yRec);
	356	hitForRec->SetZ(zRec);
30d3ea8a	357	detElemId = hitForRec->GetDetElemId();
	358	if (detElemId) iChamber = detElemId / 100 - 1;
	359	else iChamber = AliMUONConstants::ChamberNumber(zRec);
b8dc484b	360	hitForRec->SetChamberNumber(iChamber);
	361	hitForRec->SetBendingReso2(0.0);
	362	hitForRec->SetNonBendingReso2(0.0);
	363	trackParam->SetNonBendingCoor(xRec);
	364	trackParam->SetBendingCoor(yRec);
	365	trackParam->SetZ(zRec);
	366	trackParam->SetNonBendingSlope(nonBendingSlope);
	367	trackParam->SetBendingSlope(bendingSlope);
de876dd1	368	if (TMath::Abs(bendingMomentum) > 0)
48217459	369	trackParam->SetInverseBendingMomentum(1./bendingMomentum);
	370
	371	trackNew.AddHitForRecAtHit(hitForRec);
	372	trackNew.AddTrackParamAtHit(trackParam,0);
b8dc484b	373
	374	iHit1++;
	375	} // end iHit1
48217459	376
	377	trackNew.SetTrackID(track->GetTrackID());
	378	trackNew.SetTrackParamAtVertex(trackParamAtVertex);
	379	newMuonTrackRef->Add(trackNew);
b8dc484b	380
b8dc484b	381	} // end trackRef
b8dc484b	382
b8dc484b	383	delete hitForRec;
b8dc484b	384	delete trackParam;
48217459	385	return newMuonTrackRef;
b8dc484b	386	}
	387
	388	//_____________________________________________________________________________
48217459	389	AliMUONVTrackStore*
48217459	390	AliMUONRecoCheck::MakeReconstructibleTracks(const AliMUONVTrackStore& trackRefs)
b8dc484b	391	{
48217459	392	/// Calculate the number of reconstructible tracks
	393
	394	AliMUONVTrackStore* reconstructibleStore = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
b8dc484b	395
b8dc484b	396	TClonesArray* hitForRecAtHit = NULL;
	397	AliMUONHitForRec* hitForRec;
	398	Float_t zRec;
48217459	399	Int_t nTrackHits;
b8dc484b	400	Int_t isChamberInTrack[10];
	401	Int_t iChamber = 0;
	402	Bool_t isTrackOK = kTRUE;
48217459	403
	404	TIter next(trackRefs.CreateIterator());
	405	AliMUONTrack* track;
b8dc484b	406
48217459	407	while ( ( track = static_cast<AliMUONTrack*>(next()) ) )
48217459	408	{
b8dc484b	409	hitForRecAtHit = track->GetHitForRecAtHit();
	410	nTrackHits = hitForRecAtHit->GetEntriesFast();
	411	for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
48217459	412
b8dc484b	413	for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
	414	hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
	415	zRec = hitForRec->GetZ();
48217459	416	iChamber = hitForRec->GetChamberNumber();
b8dc484b	417	if (iChamber < 0 \|\| iChamber > 10) continue;
b8dc484b	418	isChamberInTrack[iChamber] = 1;
cc056de1	419	}
	420	// track is reconstructible if the particle is depositing a hit
	421	// in the following chamber combinations:
48217459	422
b8dc484b	423	isTrackOK = kTRUE;
cc056de1	424	if (!isChamberInTrack[0] && !isChamberInTrack[1]) isTrackOK = kFALSE;
	425	if (!isChamberInTrack[2] && !isChamberInTrack[3]) isTrackOK = kFALSE;
	426	if (!isChamberInTrack[4] && !isChamberInTrack[5]) isTrackOK = kFALSE;
	427	Int_t nHitsInLastStations=0;
	428	for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++)
	429	if (isChamberInTrack[ch]) nHitsInLastStations++;
	430	if(nHitsInLastStations < 3) isTrackOK = kFALSE;
48217459	431
	432	if (isTrackOK)
	433	{
	434	reconstructibleStore->Add(*track);
	435	}
b8dc484b	436	}
48217459	437
48217459	438	return reconstructibleStore;
b8dc484b	439	}
b8dc484b	440