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