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

// 13 Nov 2007:
// Added a method to create a list of reconstructed AliMUONTrack objects from
// ESD data. This is necessary since the track objects that are actually created
// during offline reconstruction are no longer stored to disk.
//  - Artur Szostak <artursz@iafrica.com>

#include "AliMUONRecoCheck.h"
#include "AliMUONRawClusterV2.h"
#include "AliMUONTrack.h"
#include "AliMUONConstants.h"
#include "AliMUONDataInterface.h"
#include "AliMUONTrackStoreV1.h"

#include "AliMCEventHandler.h"
#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliTrackReference.h"
#include "AliLog.h"
#include "AliESDEvent.h"
#include "AliESDMuonTrack.h"

#include <TParticle.h>
#include <TParticlePDG.h>
#include <Riostream.h>

#include <cassert>

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

//_____________________________________________________________________________
AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader, Char_t *pathSim)
: TObject(),
fMCEventHandler(new AliMCEventHandler()),
fDataInterface(new AliMUONDataInterface(chLoader)),
fCurrentEvent(0),
fTrackRefStore(0x0),
fRecoTrackRefStore(0x0),
fRecoTrackStore(0x0)
{
  /// Normal ctor
  fMCEventHandler->SetInputPath(pathSim);
  fMCEventHandler->InitIO("");
}

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

//_____________________________________________________________________________
void AliMUONRecoCheck::ResetStores()
{
  /// Deletes all the store objects that have been created and resets the pointers to 0x0
  delete fTrackRefStore;      fTrackRefStore = 0x0;
  delete fRecoTrackRefStore;  fRecoTrackRefStore = 0x0;
  delete fRecoTrackStore;     fRecoTrackStore = 0x0;
}

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

//_____________________________________________________________________________
AliMUONVTrackStore* AliMUONRecoCheck::ReconstructedTracks(Int_t event)
{
  /// Return the reconstructed track store for a given event
  if (!fDataInterface->IsValid()) return new AliMUONTrackStoreV1();
  
  if (event != fCurrentEvent) {
    ResetStores();
    fCurrentEvent = event;
  }
  
  if (fRecoTrackStore != 0x0) return fRecoTrackStore;
  
  fRecoTrackStore = new AliMUONTrackStoreV1();
  
  return fRecoTrackStore;
}

//_____________________________________________________________________________
AliMUONVTrackStore*
AliMUONRecoCheck::ReconstructedTracks(AliESDEvent* esd, Bool_t padMissing)
{
  /// Create and return a list of reconstructed tracks from ESD data.
  /// Note: the returned object must be deleted by the caller with the delete operator!
  /// For example:
  ///
  ///   AliMUONVTrackStore* tracks = AliMUONRecoCheck::ReconstructedTracks(esd);
  ///   delete tracks;
  ///
  /// If the padMissing flag is set to kTRUE then the missing hits on chambers [2..14]
  /// are added with the coordinate (0, 0, 0). This should be fixed in the MUON ESD
  /// structure, hopefully soon.
  
  assert( esd != NULL );
  
  AliMUONVTrackStore* store = new AliMUONTrackStoreV1;
  
  Int_t nTracks = Int_t(esd->GetNumberOfMuonTracks());
  for (Int_t i = 0; i < nTracks; i++)
  {
    AliESDMuonTrack* esdTrack = esd->GetMuonTrack(i);
    AliMUONTrack track;
    
    // Fill the track parameters at the vertex.
    AliMUONTrackParam vertexParams, ch1Params;
    vertexParams.GetParamFrom(*esdTrack);
    ch1Params.GetParamFromUncorrected(*esdTrack);
    ch1Params.GetCovFrom(*esdTrack);
    
    // Create the hit on chamber 1 (assume the same values as uncorrected track params)
    AliMUONRawClusterV2 hit(0, 0, 0); // ChamberID = 0, but unknown DE or cluster index so set those to zero.
    hit.SetXYZ(ch1Params.GetNonBendingCoor(), ch1Params.GetBendingCoor(), ch1Params.GetZ());
    hit.SetErrXY(AliMUONConstants::DefaultNonBendingReso(), AliMUONConstants::DefaultBendingReso());

    track.SetTrackParamAtVertex(&vertexParams);
    track.AddTrackParamAtCluster(ch1Params, hit, kTRUE);
    
    if (padMissing)
    {
      // We need to add fake hits so that AliMUONTrack::CompatibleTrack will be able to
      // match these tracks to the reference tracks returned by the TrackRefs method.
      for (Int_t chamber = 1; chamber < AliMUONConstants::NTrackingCh(); chamber++)
      {
        if (esdTrack->IsInMuonClusterMap(chamber))
        {
          AliMUONTrackParam fakeParams;
          fakeParams.SetZ(AliMUONConstants::DefaultChamberZ(chamber));
          fakeParams.SetBendingCoor(0);
          fakeParams.SetNonBendingCoor(0);
          AliMUONRawClusterV2 fakeHit(chamber, 0, 0); // Unknown DE or cluster index so set to zero.
          fakeHit.SetXYZ(0, 0, AliMUONConstants::DefaultChamberZ(chamber));
          fakeHit.SetErrXY(AliMUONConstants::DefaultNonBendingReso(), AliMUONConstants::DefaultBendingReso());
          track.AddTrackParamAtCluster(fakeParams, fakeHit, kTRUE);
        }
      }
      for (Int_t chamber = AliMUONConstants::NTrackingCh(); chamber < AliMUONConstants::NCh(); chamber++)
      {
        Int_t bit = 3 - (chamber - AliMUONConstants::NTrackingCh());
        if ( ((esdTrack->GetHitsPatternInTrigCh() >> bit) & 0x11) != 0 )
        {
          AliMUONTrackParam fakeParams;
          fakeParams.SetZ(AliMUONConstants::DefaultChamberZ(chamber));
          fakeParams.SetBendingCoor(0);
          fakeParams.SetNonBendingCoor(0);
          AliMUONRawClusterV2 fakeHit(chamber, 0, 0); // Unknown DE or cluster index so set to zero.
          fakeHit.SetXYZ(0, 0, AliMUONConstants::DefaultChamberZ(chamber));
          fakeHit.SetErrXY(1, 1);
          track.AddTrackParamAtCluster(fakeParams, fakeHit, kTRUE);
        }
      }
    }
    store->Add(track);
  }
  
  return store;
}

//_____________________________________________________________________________
AliMUONVTrackStore* AliMUONRecoCheck::TrackRefs(Int_t event)
{
  /// Return a track store containing the track references (converted into 
  /// MUONTrack objects) for a given event
  if (event != fCurrentEvent) {
    ResetStores();
    fCurrentEvent = event;
  }
  
  if (fTrackRefStore != 0x0) return fTrackRefStore;
  else {
    if (!fMCEventHandler->GetEvent(event)) return 0x0;
    MakeTrackRefs();
    return fTrackRefStore;
  }
}

//_____________________________________________________________________________
AliMUONVTrackStore* AliMUONRecoCheck::ReconstructibleTracks(Int_t event)
{
  /// Return a track store containing the reconstructible tracks for a given event
  if (event != fCurrentEvent) {
    ResetStores();
    fCurrentEvent = event;
  }
  
  if (fRecoTrackRefStore != 0x0) return fRecoTrackRefStore;
  else {
    if (TrackRefs(event) == 0x0) return 0x0;
    MakeReconstructibleTracks();
    return fRecoTrackRefStore;
  }
}

//_____________________________________________________________________________
void AliMUONRecoCheck::MakeTrackRefs()
{
  /// Make reconstructible tracks
  AliMUONVTrackStore *tmpTrackRefStore = new AliMUONTrackStoreV1();
  
  Double_t x, y, z, pX, pY, pZ, bendingSlope, nonBendingSlope, inverseBendingMomentum;
  TParticle* particle;
  TClonesArray* trackRefs;
  Int_t nTrackRef = fMCEventHandler->MCEvent()->GetNumberOfTracks();
  
  // loop over simulated tracks
  for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; ++iTrackRef) {
    Int_t nHits = fMCEventHandler->GetParticleAndTR(iTrackRef, particle, trackRefs);
    
    // skip empty trackRefs
    if (nHits < 1) continue;
    
    // skip trackRefs not in MUON
    AliTrackReference* trackReference0 = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(0));
    if (trackReference0->DetectorId() != AliTrackReference::kMUON) continue;
    
    AliMUONTrack track;
    
    // get track parameters at particle's vertex
    x = particle->Vx();
    y = particle->Vy();
    z = particle->Vz();
    pX = particle->Px();
    pY = particle->Py();
    pZ = particle->Pz();
    
    // compute rest of track parameters at particle's vertex
    bendingSlope = 0;
    nonBendingSlope = 0;
    inverseBendingMomentum = 0;
    if (TMath::Abs(pZ) > 0) {
      bendingSlope = pY/pZ;
      nonBendingSlope = pX/pZ;
    }
    Double_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;
    
    // set track parameters at particle's vertex
    AliMUONTrackParam trackParamAtVertex;
    trackParamAtVertex.SetNonBendingCoor(x);
    trackParamAtVertex.SetBendingCoor(y);
    trackParamAtVertex.SetZ(z);
    trackParamAtVertex.SetBendingSlope(bendingSlope);
    trackParamAtVertex.SetNonBendingSlope(nonBendingSlope);
    trackParamAtVertex.SetInverseBendingMomentum(inverseBendingMomentum);
        
    // add track parameters at vertex
    track.SetTrackParamAtVertex(&trackParamAtVertex);
    
    // loop over simulated track hits
    for (Int_t iHit = 0; iHit < nHits; ++iHit) {        
      AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(iHit));
      
      // Get track parameters of current hit
      x = trackReference->X();
      y = trackReference->Y();
      z = trackReference->Z();
      pX = trackReference->Px();
      pY = trackReference->Py();
      pZ = trackReference->Pz();
      
      // check chamberId of current trackReference
      Int_t chamberId = AliMUONConstants::ChamberNumber(z);
      if (chamberId < 0 || chamberId >= AliMUONConstants::NTrackingCh()) continue;
      
      // set hit parameters
      AliMUONRawClusterV2 hit(chamberId, 0, 0);
      hit.SetXYZ(x,y,z);
      hit.SetErrXY(0.,0.);
      
      // compute track parameters at hit
      Double_t bendingSlope = 0;
      Double_t nonBendingSlope = 0;
      Double_t inverseBendingMomentum = 0;
      if (TMath::Abs(pZ) > 0) {
	bendingSlope = pY/pZ;
	nonBendingSlope = pX/pZ;
      }
      Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
      if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
      inverseBendingMomentum *= charge;
      
      // set track parameters at hit
      AliMUONTrackParam trackParam;
      trackParam.SetNonBendingCoor(x);
      trackParam.SetBendingCoor(y);
      trackParam.SetZ(z);
      trackParam.SetBendingSlope(bendingSlope);
      trackParam.SetNonBendingSlope(nonBendingSlope);
      trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
      
      // add track parameters at current hit to the track
      track.AddTrackParamAtCluster(trackParam,hit,kTRUE);
    }
    
    track.GetTrackParamAtCluster()->Sort();
    track.SetTrackID(iTrackRef);
    tmpTrackRefStore->Add(track);
  }
  
  CleanMuonTrackRef(tmpTrackRefStore);
  
  delete tmpTrackRefStore;
}

//_____________________________________________________________________________
void AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore *tmpTrackRefStore)
{
  /// 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) 
  fTrackRefStore = new AliMUONTrackStoreV1();
  
  Double_t maxGasGap = 1.; // cm 
  Double_t x, y, z, pX, pY, pZ, x1, y1, z1, pX1, pY1, pZ1, z2;
  Double_t bendingSlope,nonBendingSlope,inverseBendingMomentum;
  
  // create iterator
  TIter next(tmpTrackRefStore->CreateIterator());
  
  // loop over tmpTrackRef
  AliMUONTrack* track;
  while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
    
    AliMUONTrack newTrack;
    
    // loop over tmpTrackRef's hits
    Int_t iHit1 = 0;
    Int_t nTrackHits = track->GetNClusters();
    while (iHit1 < nTrackHits) {
      AliMUONTrackParam *trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit1);
      
      // get track parameters at hit1
      x1  = trackParam1->GetNonBendingCoor();
      y1  = trackParam1->GetBendingCoor();
      z1  = trackParam1->GetZ();
      pX1 = trackParam1->Px();
      pY1 = trackParam1->Py();
      pZ1 = trackParam1->Pz();
      
      // prepare new track parameters
      x  = x1;
      y  = y1;
      z  = z1;
      pX = pX1;
      pY = pY1;
      pZ = pZ1;
      
      // loop over next tmpTrackRef's hits
      Int_t nCombinedHits = 1;
      for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
	AliMUONTrackParam *trackParam2 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit2);
        
	// get z position of hit2
	z2 = trackParam2->GetZ();
	
	// complete new track parameters if hit2 is on the same detection element
        if ( TMath::Abs(z2-z1) < maxGasGap ) {
	  x  += trackParam2->GetNonBendingCoor();
	  y  += trackParam2->GetBendingCoor();
	  z  += z2;
	  pX += trackParam2->Px();
	  pY += trackParam2->Py();
	  pZ += trackParam2->Pz();
	  nCombinedHits++;
          iHit1 = iHit2;
        }
        
      }
      
      // finalize new track parameters
      x  /= (Double_t)nCombinedHits;
      y  /= (Double_t)nCombinedHits;
      z  /= (Double_t)nCombinedHits;
      pX /= (Double_t)nCombinedHits;
      pY /= (Double_t)nCombinedHits;
      pZ /= (Double_t)nCombinedHits;
      bendingSlope = 0;
      nonBendingSlope = 0;
      inverseBendingMomentum = 0;
      if (TMath::Abs(pZ) > 0) {
	bendingSlope = pY/pZ;
	nonBendingSlope = pX/pZ;
      }
      Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
      if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
      inverseBendingMomentum *= trackParam1->GetCharge();
      
      // set hit parameters
      AliMUONRawClusterV2 hit(trackParam1->GetClusterPtr()->GetChamberId(), 0, 0);
      hit.SetXYZ(x,y,z);
      hit.SetErrXY(0.,0.);
      
      // set new track parameters at new hit
      AliMUONTrackParam trackParam;
      trackParam.SetNonBendingCoor(x);
      trackParam.SetBendingCoor(y);
      trackParam.SetZ(z);
      trackParam.SetBendingSlope(bendingSlope);
      trackParam.SetNonBendingSlope(nonBendingSlope);
      trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
      
      // add track parameters at current hit to the track
      newTrack.AddTrackParamAtCluster(trackParam,hit,kTRUE);
      
      iHit1++;
    }
    
    newTrack.GetTrackParamAtCluster()->Sort();
    newTrack.SetTrackID(track->GetTrackID());
    newTrack.SetTrackParamAtVertex(track->GetTrackParamAtVertex());
    fTrackRefStore->Add(newTrack);
    
  }
  
}

//_____________________________________________________________________________
void AliMUONRecoCheck::MakeReconstructibleTracks()
{
  /// Calculate the number of reconstructible tracks
  fRecoTrackRefStore = new AliMUONTrackStoreV1();
  
  // create iterator on trackRef
  TIter next(fTrackRefStore->CreateIterator());
  
  // loop over trackRef
  AliMUONTrack* track;
  while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
    
    Bool_t* chamberInTrack = new Bool_t(AliMUONConstants::NTrackingCh());
    for (Int_t iCh = 0; iCh < AliMUONConstants::NTrackingCh(); iCh++) chamberInTrack[iCh] = kFALSE;
    
    // loop over trackRef's hits to get hit chambers
    Int_t nTrackHits = track->GetNClusters();
    for (Int_t iHit = 0; iHit < nTrackHits; iHit++) {
      AliMUONVCluster* hit = ((AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit))->GetClusterPtr(); 
      chamberInTrack[hit->GetChamberId()] = kTRUE;
    } 
    
    // track is reconstructible if the particle is depositing a hit
    // in the following chamber combinations:
    Bool_t trackOK = kTRUE;
    if (!chamberInTrack[0] && !chamberInTrack[1]) trackOK = kFALSE;
    if (!chamberInTrack[2] && !chamberInTrack[3]) trackOK = kFALSE;
    if (!chamberInTrack[4] && !chamberInTrack[5]) trackOK = kFALSE;
    Int_t nHitsInLastStations = 0;
    for (Int_t iCh = 6; iCh < AliMUONConstants::NTrackingCh(); iCh++)
      if (chamberInTrack[iCh]) nHitsInLastStations++; 
    if(nHitsInLastStations < 3) trackOK = kFALSE;
    
    // Add reconstructible tracks to fRecoTrackRefStore
    if (trackOK) fRecoTrackRefStore->Add(*track);
    
    delete [] chamberInTrack;
  }

}
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
9bdda5f6	26	// 13 Nov 2007:
	27	// Added a method to create a list of reconstructed AliMUONTrack objects from
	28	// ESD data. This is necessary since the track objects that are actually created
	29	// during offline reconstruction are no longer stored to disk.
	30	// - Artur Szostak <artursz@iafrica.com>
	31
cc056de1	32	#include "AliMUONRecoCheck.h"
96ebe67e	33	#include "AliMUONRawClusterV2.h"
cc056de1	34	#include "AliMUONTrack.h"
cc056de1	35	#include "AliMUONConstants.h"
48217459	36	#include "AliMUONDataInterface.h"
96ebe67e	37	#include "AliMUONTrackStoreV1.h"
	38
	39	#include "AliMCEventHandler.h"
	40	#include "AliMCEvent.h"
b8dc484b	41	#include "AliStack.h"
cc056de1	42	#include "AliTrackReference.h"
9bdda5f6	43	#include "AliLog.h"
	44	#include "AliESDEvent.h"
	45	#include "AliESDMuonTrack.h"
96ebe67e	46
cc056de1	47	#include <TParticle.h>
5b1afea0	48	#include <TParticlePDG.h>
96ebe67e	49	#include <Riostream.h>
96ebe67e	50
9bdda5f6	51	#include <cassert>
9bdda5f6	52
5398f946	53	/// \cond CLASSIMP
b8dc484b	54	ClassImp(AliMUONRecoCheck)
5398f946	55	/// \endcond
b8dc484b	56
068e4c36	57	//_____________________________________________________________________________
96ebe67e	58	AliMUONRecoCheck::AliMUONRecoCheck(Char_t chLoader, Char_t pathSim)
48217459	59	: TObject(),
96ebe67e	60	fMCEventHandler(new AliMCEventHandler()),
	61	fDataInterface(new AliMUONDataInterface(chLoader)),
	62	fCurrentEvent(0),
	63	fTrackRefStore(0x0),
	64	fRecoTrackRefStore(0x0),
	65	fRecoTrackStore(0x0)
068e4c36	66	{
48217459	67	/// Normal ctor
96ebe67e	68	fMCEventHandler->SetInputPath(pathSim);
96ebe67e	69	fMCEventHandler->InitIO("");
068e4c36	70	}
068e4c36	71
b8dc484b	72	//_____________________________________________________________________________
48217459	73	AliMUONRecoCheck::~AliMUONRecoCheck()
b8dc484b	74	{
48217459	75	/// Destructor
96ebe67e	76	delete fMCEventHandler;
48217459	77	delete fDataInterface;
96ebe67e	78	ResetStores();
	79	}
	80
	81	//_____________________________________________________________________________
	82	void AliMUONRecoCheck::ResetStores()
	83	{
	84	/// Deletes all the store objects that have been created and resets the pointers to 0x0
	85	delete fTrackRefStore; fTrackRefStore = 0x0;
	86	delete fRecoTrackRefStore; fRecoTrackRefStore = 0x0;
	87	delete fRecoTrackStore; fRecoTrackStore = 0x0;
	88	}
	89
	90	//_____________________________________________________________________________
	91	Int_t AliMUONRecoCheck::NumberOfEvents() const
	92	{
	93	/// Return the number of events
	94	if (fDataInterface->IsValid()) return fDataInterface->NumberOfEvents();
	95	return 0;
b8dc484b	96	}
b8dc484b	97
48217459	98	//_____________________________________________________________________________
96ebe67e	99	AliMUONVTrackStore* AliMUONRecoCheck::ReconstructedTracks(Int_t event)
48217459	100	{
48217459	101	/// Return the reconstructed track store for a given event
96ebe67e	102	if (!fDataInterface->IsValid()) return new AliMUONTrackStoreV1();
	103
	104	if (event != fCurrentEvent) {
	105	ResetStores();
	106	fCurrentEvent = event;
	107	}
	108
	109	if (fRecoTrackStore != 0x0) return fRecoTrackStore;
	110
	111	fRecoTrackStore = new AliMUONTrackStoreV1();
	112
	113	return fRecoTrackStore;
48217459	114	}
45a05f07	115
b8dc484b	116	//_____________________________________________________________________________
9bdda5f6	117	AliMUONVTrackStore*
	118	AliMUONRecoCheck::ReconstructedTracks(AliESDEvent* esd, Bool_t padMissing)
	119	{
	120	/// Create and return a list of reconstructed tracks from ESD data.
	121	/// Note: the returned object must be deleted by the caller with the delete operator!
	122	/// For example:
	123	///
	124	/// AliMUONVTrackStore* tracks = AliMUONRecoCheck::ReconstructedTracks(esd);
	125	/// delete tracks;
	126	///
	127	/// If the padMissing flag is set to kTRUE then the missing hits on chambers [2..14]
	128	/// are added with the coordinate (0, 0, 0). This should be fixed in the MUON ESD
	129	/// structure, hopefully soon.
	130
	131	assert( esd != NULL );
	132
	133	AliMUONVTrackStore* store = new AliMUONTrackStoreV1;
	134
	135	Int_t nTracks = Int_t(esd->GetNumberOfMuonTracks());
	136	for (Int_t i = 0; i < nTracks; i++)
	137	{
	138	AliESDMuonTrack* esdTrack = esd->GetMuonTrack(i);
	139	AliMUONTrack track;
	140
	141	// Fill the track parameters at the vertex.
	142	AliMUONTrackParam vertexParams, ch1Params;
	143	vertexParams.GetParamFrom(*esdTrack);
	144	ch1Params.GetParamFromUncorrected(*esdTrack);
	145	ch1Params.GetCovFrom(*esdTrack);
	146
	147	// Create the hit on chamber 1 (assume the same values as uncorrected track params)
	148	AliMUONRawClusterV2 hit(0, 0, 0); // ChamberID = 0, but unknown DE or cluster index so set those to zero.
	149	hit.SetXYZ(ch1Params.GetNonBendingCoor(), ch1Params.GetBendingCoor(), ch1Params.GetZ());
	150	hit.SetErrXY(AliMUONConstants::DefaultNonBendingReso(), AliMUONConstants::DefaultBendingReso());
	151
	152	track.SetTrackParamAtVertex(&vertexParams);
	153	track.AddTrackParamAtCluster(ch1Params, hit, kTRUE);
	154
	155	if (padMissing)
	156	{
	157	// We need to add fake hits so that AliMUONTrack::CompatibleTrack will be able to
	158	// match these tracks to the reference tracks returned by the TrackRefs method.
	159	for (Int_t chamber = 1; chamber < AliMUONConstants::NTrackingCh(); chamber++)
	160	{
	161	if (esdTrack->IsInMuonClusterMap(chamber))
	162	{
	163	AliMUONTrackParam fakeParams;
	164	fakeParams.SetZ(AliMUONConstants::DefaultChamberZ(chamber));
	165	fakeParams.SetBendingCoor(0);
	166	fakeParams.SetNonBendingCoor(0);
	167	AliMUONRawClusterV2 fakeHit(chamber, 0, 0); // Unknown DE or cluster index so set to zero.
	168	fakeHit.SetXYZ(0, 0, AliMUONConstants::DefaultChamberZ(chamber));
	169	fakeHit.SetErrXY(AliMUONConstants::DefaultNonBendingReso(), AliMUONConstants::DefaultBendingReso());
	170	track.AddTrackParamAtCluster(fakeParams, fakeHit, kTRUE);
	171	}
	172	}
	173	for (Int_t chamber = AliMUONConstants::NTrackingCh(); chamber < AliMUONConstants::NCh(); chamber++)
	174	{
	175	Int_t bit = 3 - (chamber - AliMUONConstants::NTrackingCh());
	176	if ( ((esdTrack->GetHitsPatternInTrigCh() >> bit) & 0x11) != 0 )
	177	{
	178	AliMUONTrackParam fakeParams;
	179	fakeParams.SetZ(AliMUONConstants::DefaultChamberZ(chamber));
	180	fakeParams.SetBendingCoor(0);
181	fakeParams.SetNonBendingCoor(0);
182	AliMUONRawClusterV2 fakeHit(chamber, 0, 0); // Unknown DE or cluster index so set to zero.
183	fakeHit.SetXYZ(0, 0, AliMUONConstants::DefaultChamberZ(chamber));
184	fakeHit.SetErrXY(1, 1);
185	track.AddTrackParamAtCluster(fakeParams, fakeHit, kTRUE);
186	}
187	}
188	}
189	store->Add(track);
190	}
191
192	return store;
193	}
194
195	//_____________________________________________________________________________
96ebe67e	196	AliMUONVTrackStore* AliMUONRecoCheck::TrackRefs(Int_t event)
b8dc484b	197	{
48217459	198	/// Return a track store containing the track references (converted into
48217459	199	/// MUONTrack objects) for a given event
96ebe67e	200	if (event != fCurrentEvent) {
	201	ResetStores();
	202	fCurrentEvent = event;
	203	}
	204
	205	if (fTrackRefStore != 0x0) return fTrackRefStore;
	206	else {
	207	if (!fMCEventHandler->GetEvent(event)) return 0x0;
	208	MakeTrackRefs();
	209	return fTrackRefStore;
	210	}
b8dc484b	211	}
	212
	213	//_____________________________________________________________________________
96ebe67e	214	AliMUONVTrackStore* AliMUONRecoCheck::ReconstructibleTracks(Int_t event)
b8dc484b	215	{
48217459	216	/// Return a track store containing the reconstructible tracks for a given event
96ebe67e	217	if (event != fCurrentEvent) {
	218	ResetStores();
	219	fCurrentEvent = event;
	220	}
	221
	222	if (fRecoTrackRefStore != 0x0) return fRecoTrackRefStore;
	223	else {
	224	if (TrackRefs(event) == 0x0) return 0x0;
	225	MakeReconstructibleTracks();
	226	return fRecoTrackRefStore;
	227	}
48217459	228	}
b8dc484b	229
48217459	230	//_____________________________________________________________________________
96ebe67e	231	void AliMUONRecoCheck::MakeTrackRefs()
48217459	232	{
48217459	233	/// Make reconstructible tracks
96ebe67e	234	AliMUONVTrackStore *tmpTrackRefStore = new AliMUONTrackStoreV1();
48217459	235
96ebe67e	236	Double_t x, y, z, pX, pY, pZ, bendingSlope, nonBendingSlope, inverseBendingMomentum;
	237	TParticle* particle;
	238	TClonesArray* trackRefs;
	239	Int_t nTrackRef = fMCEventHandler->MCEvent()->GetNumberOfTracks();
b8dc484b	240
96ebe67e	241	// loop over simulated tracks
	242	for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; ++iTrackRef) {
	243	Int_t nHits = fMCEventHandler->GetParticleAndTR(iTrackRef, particle, trackRefs);
48217459	244
96ebe67e	245	// skip empty trackRefs
	246	if (nHits < 1) continue;
	247
	248	// skip trackRefs not in MUON
	249	AliTrackReference* trackReference0 = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(0));
	250	if (trackReference0->DetectorId() != AliTrackReference::kMUON) continue;
	251
	252	AliMUONTrack track;
b8dc484b	253
96ebe67e	254	// get track parameters at particle's vertex
	255	x = particle->Vx();
	256	y = particle->Vy();
	257	z = particle->Vz();
	258	pX = particle->Px();
	259	pY = particle->Py();
	260	pZ = particle->Pz();
dba3c809	261
96ebe67e	262	// compute rest of track parameters at particle's vertex
	263	bendingSlope = 0;
	264	nonBendingSlope = 0;
	265	inverseBendingMomentum = 0;
	266	if (TMath::Abs(pZ) > 0) {
	267	bendingSlope = pY/pZ;
	268	nonBendingSlope = pX/pZ;
	269	}
	270	Double_t pYZ = TMath::Sqrt(pYpY+pZpZ);
	271	if (pYZ >0) inverseBendingMomentum = 1/pYZ;
	272	TParticlePDG* ppdg = particle->GetPDG(1);
	273	Int_t charge = (Int_t)(ppdg->Charge()/3.0);
	274	inverseBendingMomentum *= charge;
	275
	276	// set track parameters at particle's vertex
	277	AliMUONTrackParam trackParamAtVertex;
	278	trackParamAtVertex.SetNonBendingCoor(x);
	279	trackParamAtVertex.SetBendingCoor(y);
	280	trackParamAtVertex.SetZ(z);
	281	trackParamAtVertex.SetBendingSlope(bendingSlope);
	282	trackParamAtVertex.SetNonBendingSlope(nonBendingSlope);
	283	trackParamAtVertex.SetInverseBendingMomentum(inverseBendingMomentum);
	284
	285	// add track parameters at vertex
	286	track.SetTrackParamAtVertex(&trackParamAtVertex);
48217459	287
96ebe67e	288	// loop over simulated track hits
	289	for (Int_t iHit = 0; iHit < nHits; ++iHit) {
	290	AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(iHit));
068e4c36	291
96ebe67e	292	// Get track parameters of current hit
	293	x = trackReference->X();
	294	y = trackReference->Y();
	295	z = trackReference->Z();
	296	pX = trackReference->Px();
	297	pY = trackReference->Py();
	298	pZ = trackReference->Pz();
b8dc484b	299
96ebe67e	300	// check chamberId of current trackReference
	301	Int_t chamberId = AliMUONConstants::ChamberNumber(z);
	302	if (chamberId < 0 \|\| chamberId >= AliMUONConstants::NTrackingCh()) continue;
48217459	303
96ebe67e	304	// set hit parameters
	305	AliMUONRawClusterV2 hit(chamberId, 0, 0);
	306	hit.SetXYZ(x,y,z);
	307	hit.SetErrXY(0.,0.);
b8dc484b	308
96ebe67e	309	// compute track parameters at hit
	310	Double_t bendingSlope = 0;
	311	Double_t nonBendingSlope = 0;
	312	Double_t inverseBendingMomentum = 0;
	313	if (TMath::Abs(pZ) > 0) {
	314	bendingSlope = pY/pZ;
	315	nonBendingSlope = pX/pZ;
b8dc484b	316	}
96ebe67e	317	Double_t pYZ = TMath::Sqrt(pYpY+pZpZ);
	318	if (pYZ >0) inverseBendingMomentum = 1/pYZ;
	319	inverseBendingMomentum *= charge;
	320
	321	// set track parameters at hit
	322	AliMUONTrackParam trackParam;
	323	trackParam.SetNonBendingCoor(x);
	324	trackParam.SetBendingCoor(y);
	325	trackParam.SetZ(z);
	326	trackParam.SetBendingSlope(bendingSlope);
	327	trackParam.SetNonBendingSlope(nonBendingSlope);
	328	trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
	329
	330	// add track parameters at current hit to the track
	331	track.AddTrackParamAtCluster(trackParam,hit,kTRUE);
	332	}
	333
	334	track.GetTrackParamAtCluster()->Sort();
	335	track.SetTrackID(iTrackRef);
	336	tmpTrackRefStore->Add(track);
b8dc484b	337	}
b8dc484b	338
96ebe67e	339	CleanMuonTrackRef(tmpTrackRefStore);
48217459	340
96ebe67e	341	delete tmpTrackRefStore;
b8dc484b	342	}
	343
	344	//_____________________________________________________________________________
96ebe67e	345	void AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore *tmpTrackRefStore)
b8dc484b	346	{
48217459	347	/// Re-calculate hits parameters because two AliTrackReferences are recorded for
	348	/// each chamber (one when particle is entering + one when particle is leaving
	349	/// the sensitive volume)
96ebe67e	350	fTrackRefStore = new AliMUONTrackStoreV1();
48217459	351
96ebe67e	352	Double_t maxGasGap = 1.; // cm
	353	Double_t x, y, z, pX, pY, pZ, x1, y1, z1, pX1, pY1, pZ1, z2;
	354	Double_t bendingSlope,nonBendingSlope,inverseBendingMomentum;
48217459	355
96ebe67e	356	// create iterator
96ebe67e	357	TIter next(tmpTrackRefStore->CreateIterator());
48217459	358
96ebe67e	359	// loop over tmpTrackRef
48217459	360	AliMUONTrack* track;
96ebe67e	361	while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
	362
	363	AliMUONTrack newTrack;
	364
	365	// loop over tmpTrackRef's hits
	366	Int_t iHit1 = 0;
	367	Int_t nTrackHits = track->GetNClusters();
	368	while (iHit1 < nTrackHits) {
	369	AliMUONTrackParam trackParam1 = (AliMUONTrackParam) track->GetTrackParamAtCluster()->UncheckedAt(iHit1);
	370
	371	// get track parameters at hit1
	372	x1 = trackParam1->GetNonBendingCoor();
	373	y1 = trackParam1->GetBendingCoor();
	374	z1 = trackParam1->GetZ();
	375	pX1 = trackParam1->Px();
	376	pY1 = trackParam1->Py();
	377	pZ1 = trackParam1->Pz();
	378
	379	// prepare new track parameters
	380	x = x1;
	381	y = y1;
	382	z = z1;
	383	pX = pX1;
	384	pY = pY1;
	385	pZ = pZ1;
	386
	387	// loop over next tmpTrackRef's hits
	388	Int_t nCombinedHits = 1;
	389	for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
	390	AliMUONTrackParam trackParam2 = (AliMUONTrackParam) track->GetTrackParamAtCluster()->UncheckedAt(iHit2);
48217459	391
96ebe67e	392	// get z position of hit2
	393	z2 = trackParam2->GetZ();
	394
	395	// complete new track parameters if hit2 is on the same detection element
	396	if ( TMath::Abs(z2-z1) < maxGasGap ) {
	397	x += trackParam2->GetNonBendingCoor();
	398	y += trackParam2->GetBendingCoor();
	399	z += z2;
	400	pX += trackParam2->Px();
	401	pY += trackParam2->Py();
	402	pZ += trackParam2->Pz();
	403	nCombinedHits++;
48217459	404	iHit1 = iHit2;
	405	}
	406
96ebe67e	407	}
	408
	409	// finalize new track parameters
	410	x /= (Double_t)nCombinedHits;
	411	y /= (Double_t)nCombinedHits;
	412	z /= (Double_t)nCombinedHits;
	413	pX /= (Double_t)nCombinedHits;
	414	pY /= (Double_t)nCombinedHits;
	415	pZ /= (Double_t)nCombinedHits;
	416	bendingSlope = 0;
	417	nonBendingSlope = 0;
	418	inverseBendingMomentum = 0;
	419	if (TMath::Abs(pZ) > 0) {
	420	bendingSlope = pY/pZ;
	421	nonBendingSlope = pX/pZ;
	422	}
	423	Double_t pYZ = TMath::Sqrt(pYpY+pZpZ);
	424	if (pYZ >0) inverseBendingMomentum = 1/pYZ;
	425	inverseBendingMomentum *= trackParam1->GetCharge();
b8dc484b	426
96ebe67e	427	// set hit parameters
	428	AliMUONRawClusterV2 hit(trackParam1->GetClusterPtr()->GetChamberId(), 0, 0);
	429	hit.SetXYZ(x,y,z);
	430	hit.SetErrXY(0.,0.);
48217459	431
96ebe67e	432	// set new track parameters at new hit
	433	AliMUONTrackParam trackParam;
	434	trackParam.SetNonBendingCoor(x);
	435	trackParam.SetBendingCoor(y);
	436	trackParam.SetZ(z);
	437	trackParam.SetBendingSlope(bendingSlope);
	438	trackParam.SetNonBendingSlope(nonBendingSlope);
	439	trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
	440
	441	// add track parameters at current hit to the track
	442	newTrack.AddTrackParamAtCluster(trackParam,hit,kTRUE);
b8dc484b	443
b8dc484b	444	iHit1++;
96ebe67e	445	}
48217459	446
96ebe67e	447	newTrack.GetTrackParamAtCluster()->Sort();
	448	newTrack.SetTrackID(track->GetTrackID());
	449	newTrack.SetTrackParamAtVertex(track->GetTrackParamAtVertex());
	450	fTrackRefStore->Add(newTrack);
b8dc484b	451
96ebe67e	452	}
b8dc484b	453
b8dc484b	454	}
	455
	456	//_____________________________________________________________________________
96ebe67e	457	void AliMUONRecoCheck::MakeReconstructibleTracks()
b8dc484b	458	{
48217459	459	/// Calculate the number of reconstructible tracks
96ebe67e	460	fRecoTrackRefStore = new AliMUONTrackStoreV1();
48217459	461
96ebe67e	462	// create iterator on trackRef
96ebe67e	463	TIter next(fTrackRefStore->CreateIterator());
b8dc484b	464
96ebe67e	465	// loop over trackRef
48217459	466	AliMUONTrack* track;
96ebe67e	467	while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
	468
	469	Bool_t* chamberInTrack = new Bool_t(AliMUONConstants::NTrackingCh());
	470	for (Int_t iCh = 0; iCh < AliMUONConstants::NTrackingCh(); iCh++) chamberInTrack[iCh] = kFALSE;
48217459	471
96ebe67e	472	// loop over trackRef's hits to get hit chambers
	473	Int_t nTrackHits = track->GetNClusters();
	474	for (Int_t iHit = 0; iHit < nTrackHits; iHit++) {
	475	AliMUONVCluster* hit = ((AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit))->GetClusterPtr();
	476	chamberInTrack[hit->GetChamberId()] = kTRUE;
cc056de1	477	}
96ebe67e	478
cc056de1	479	// track is reconstructible if the particle is depositing a hit
cc056de1	480	// in the following chamber combinations:
96ebe67e	481	Bool_t trackOK = kTRUE;
	482	if (!chamberInTrack[0] && !chamberInTrack[1]) trackOK = kFALSE;
	483	if (!chamberInTrack[2] && !chamberInTrack[3]) trackOK = kFALSE;
	484	if (!chamberInTrack[4] && !chamberInTrack[5]) trackOK = kFALSE;
	485	Int_t nHitsInLastStations = 0;
	486	for (Int_t iCh = 6; iCh < AliMUONConstants::NTrackingCh(); iCh++)
	487	if (chamberInTrack[iCh]) nHitsInLastStations++;
	488	if(nHitsInLastStations < 3) trackOK = kFALSE;
48217459	489
96ebe67e	490	// Add reconstructible tracks to fRecoTrackRefStore
96ebe67e	491	if (trackOK) fRecoTrackRefStore->Add(*track);
48217459	492
96ebe67e	493	delete [] chamberInTrack;
b8dc484b	494	}
48217459	495
b8dc484b	496	}
b8dc484b	497