[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>
// 25 Jan 2008:
// Use the new ESDInterface to create MUON objects from ESD data
// - Philippe Pillot
// 

#include "AliMUONRecoCheck.h"
#include "AliMUONTrack.h"
#include "AliMUONVTrackStore.h"
#include "AliMUONVCluster.h"
#include "AliMUONVClusterStore.h"
#include "AliMUONConstants.h"
#include "AliMUONESDInterface.h"
#include "AliMUONTrackParam.h"
#include "AliMCEventHandler.h"
#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliTrackReference.h"
#include "AliLog.h"
#include "AliESDEvent.h"
#include "AliESDMuonTrack.h"

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

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

//_____________________________________________________________________________
AliMUONRecoCheck::AliMUONRecoCheck(const Char_t *esdFileName, const Char_t *pathSim)
: TObject(),
fMCEventHandler(new AliMCEventHandler()),
fESDEvent(new AliESDEvent()),
fESDTree (0x0),
fESDFile (0x0),
fCurrentEvent(0),
fTrackRefStore(0x0),
fRecoTrackRefStore(0x0),
fRecoTrackStore(0x0),
fESDEventOwner(kTRUE)
{
  /// Normal ctor
  
  // TrackRefs and Particules
  fMCEventHandler->SetInputPath(pathSim);
  fMCEventHandler->InitIO("");
  
  // ESD MUON Tracks
  fESDFile = TFile::Open(esdFileName); // open the file
  if (!fESDFile || !fESDFile->IsOpen()) {
    AliError(Form("opening ESD file %s failed", esdFileName));
    fESDFile = 0x0;
    return;
  }
  fESDTree = (TTree*) fESDFile->Get("esdTree"); // get the tree
  if (!fESDTree) {
    AliError("no ESD tree found");
    fESDFile->Close();
    fESDFile = 0x0;
    return;
  }
  fESDEvent->ReadFromTree(fESDTree); // link fESDEvent to the tree
}

//_____________________________________________________________________________
AliMUONRecoCheck::AliMUONRecoCheck(AliESDEvent *esdEvent, AliMCEventHandler *mcEventHandler)
: TObject(),
fMCEventHandler(0),
fESDEvent(0),
fESDTree (0x0),
fESDFile (0x0),
fCurrentEvent(0),
fTrackRefStore(0x0),
fRecoTrackRefStore(0x0),
fRecoTrackStore(0x0),
fESDEventOwner(kFALSE)
{
  /// Normal ctor
  
  // TrackRefs and Particules
  fMCEventHandler = mcEventHandler;
  
  // ESD MUON Tracks
  fESDEvent = esdEvent;
  
}

//_____________________________________________________________________________
AliMUONRecoCheck::~AliMUONRecoCheck()
{
  /// Destructor
  if (fESDEventOwner) {
    delete fMCEventHandler;
    delete fESDEvent;
    if (fESDFile) fESDFile->Close();
  }
  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 (fESDEventOwner && fESDTree) return fESDTree->GetEntries();
  return 0;
}

//_____________________________________________________________________________
AliMUONVTrackStore* AliMUONRecoCheck::ReconstructedTracks(Int_t event)
{
  /// Return a track store containing the reconstructed tracks (converted into 
  /// MUONTrack objects) for a given event
  
  if (!fESDEventOwner) {
    MakeReconstructedTracks();
    return fRecoTrackStore;
  }

  if (event != fCurrentEvent) {
    ResetStores();
    fCurrentEvent = event;
  }
  
  if (fRecoTrackStore != 0x0) return fRecoTrackStore;
  else {
    if (!fESDTree) return 0x0;
    if (fESDTree->GetEvent(event) <= 0) {
      AliError(Form("fails to read ESD object for event %d", event));
      return 0x0;
    }
    MakeReconstructedTracks();
    return fRecoTrackStore;
  }
}

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

  if (event != fCurrentEvent) {
    ResetStores();
    fCurrentEvent = event;
  }
  
  if (fTrackRefStore != 0x0) return fTrackRefStore;
  else {
    if (!fMCEventHandler->GetEvent(event)) {
      AliError(Form("fails to read MC objects for event %d", 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 (!fESDEventOwner) {
    if (TrackRefs(event) == 0x0) return 0x0;
    MakeReconstructibleTracks();
    return fRecoTrackRefStore;
  }

  if (event != fCurrentEvent) {
    ResetStores();
    fCurrentEvent = event;
  }
  
  if (fRecoTrackRefStore != 0x0) return fRecoTrackRefStore;
  else {
    if (TrackRefs(event) == 0x0) return 0x0;
    MakeReconstructibleTracks();
    return fRecoTrackRefStore;
  }
}

//_____________________________________________________________________________
void AliMUONRecoCheck::MakeReconstructedTracks()
{
  /// Make reconstructed tracks
  if (!(fRecoTrackStore = AliMUONESDInterface::NewTrackStore())) return;
  
  // loop over all reconstructed tracks and add them to the store (skip ghosts)
  Int_t nTracks = (Int_t) fESDEvent->GetNumberOfMuonTracks();
  for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
    AliESDMuonTrack* esdTrack = fESDEvent->GetMuonTrack(iTrack);
    if (esdTrack->ContainTrackerData()) AliMUONESDInterface::Add(*esdTrack, *fRecoTrackStore);
  }
  
}

//_____________________________________________________________________________
void AliMUONRecoCheck::MakeTrackRefs()
{
  /// Make reconstructible tracks
  AliMUONVTrackStore *tmpTrackRefStore = AliMUONESDInterface::NewTrackStore();
  if (!tmpTrackRefStore) return;
  
  Double_t x, y, z, pX, pY, pZ, bendingSlope, nonBendingSlope, inverseBendingMomentum;
  TParticle* particle;
  TClonesArray* trackRefs;
  Int_t nTrackRef = fMCEventHandler->MCEvent()->GetNumberOfTracks();
  AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
  if (!cStore) return;
  AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
  
  // 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;
    
    // get the particle charge for further calculation
    TParticlePDG* ppdg = particle->GetPDG();
    Int_t charge = ppdg != NULL ? (Int_t)(ppdg->Charge()/3.0) : 0;
    
    AliMUONTrack track;
    
    // loop over simulated track hits
    for (Int_t iHit = 0; iHit < nHits; ++iHit) {        
      AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(iHit));
      
      // skip trackRefs not in MUON
      if (trackReference->DetectorId() != AliTrackReference::kMUON) continue;
    
      // 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 detElemId = trackReference->UserId();
      Int_t chamberId = detElemId / 100 - 1;
      if (chamberId < 0 || chamberId >= AliMUONConstants::NTrackingCh()) continue;
      
      // set hit parameters
      hit->SetUniqueID(AliMUONVCluster::BuildUniqueID(chamberId, detElemId, iHit));
      hit->SetXYZ(x,y,z);
      hit->SetErrXY(0.,0.);
      
      // compute track parameters at hit
      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 *= 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);
    }
    
    // if none of the track hits was in MUON, goto the next track
    if (track.GetNClusters() < 1) continue;
    
    // 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;
    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);
    
    // store the track
    track.SetUniqueID(iTrackRef);
    tmpTrackRefStore->Add(track);
  }
  
  CleanMuonTrackRef(tmpTrackRefStore);
  
  delete hit;
  delete cStore;
  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) 
  if (!(fTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
  
  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;
  AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
  if (!cStore) return;
  AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
  
  // 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
      hit->SetUniqueID(trackParam1->GetClusterPtr()->GetUniqueID());
      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.SetUniqueID(track->GetUniqueID());
    newTrack.SetTrackParamAtVertex(track->GetTrackParamAtVertex());
    fTrackRefStore->Add(newTrack);
    
  }
  
  delete hit;
  delete cStore;
}

//_____________________________________________________________________________
void AliMUONRecoCheck::MakeReconstructibleTracks()
{
  /// Isolate the reconstructible tracks
  if (!(fRecoTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
  
  // 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>
103e6575	31	// 25 Jan 2008:
	32	// Use the new ESDInterface to create MUON objects from ESD data
	33	// - Philippe Pillot
7c29acf1	34	//
9bdda5f6	35
cc056de1	36	#include "AliMUONRecoCheck.h"
cc056de1	37	#include "AliMUONTrack.h"
103e6575	38	#include "AliMUONVTrackStore.h"
	39	#include "AliMUONVCluster.h"
	40	#include "AliMUONVClusterStore.h"
cc056de1	41	#include "AliMUONConstants.h"
103e6575	42	#include "AliMUONESDInterface.h"
933daf4c	43	#include "AliMUONTrackParam.h"
96ebe67e	44	#include "AliMCEventHandler.h"
96ebe67e	45	#include "AliMCEvent.h"
b8dc484b	46	#include "AliStack.h"
cc056de1	47	#include "AliTrackReference.h"
9bdda5f6	48	#include "AliLog.h"
	49	#include "AliESDEvent.h"
	50	#include "AliESDMuonTrack.h"
96ebe67e	51
61fed964	52	#include <TFile.h>
61fed964	53	#include <TTree.h>
cc056de1	54	#include <TParticle.h>
5b1afea0	55	#include <TParticlePDG.h>
96ebe67e	56	#include <Riostream.h>
96ebe67e	57
5398f946	58	/// \cond CLASSIMP
b8dc484b	59	ClassImp(AliMUONRecoCheck)
5398f946	60	/// \endcond
b8dc484b	61
068e4c36	62	//_____________________________________________________________________________
1f2d22de	63	AliMUONRecoCheck::AliMUONRecoCheck(const Char_t esdFileName, const Char_t pathSim)
48217459	64	: TObject(),
96ebe67e	65	fMCEventHandler(new AliMCEventHandler()),
61fed964	66	fESDEvent(new AliESDEvent()),
	67	fESDTree (0x0),
	68	fESDFile (0x0),
96ebe67e	69	fCurrentEvent(0),
	70	fTrackRefStore(0x0),
	71	fRecoTrackRefStore(0x0),
7c29acf1	72	fRecoTrackStore(0x0),
7c29acf1	73	fESDEventOwner(kTRUE)
068e4c36	74	{
48217459	75	/// Normal ctor
61fed964	76
61fed964	77	// TrackRefs and Particules
96ebe67e	78	fMCEventHandler->SetInputPath(pathSim);
96ebe67e	79	fMCEventHandler->InitIO("");
61fed964	80
	81	// ESD MUON Tracks
	82	fESDFile = TFile::Open(esdFileName); // open the file
	83	if (!fESDFile \|\| !fESDFile->IsOpen()) {
	84	AliError(Form("opening ESD file %s failed", esdFileName));
	85	fESDFile = 0x0;
	86	return;
	87	}
	88	fESDTree = (TTree*) fESDFile->Get("esdTree"); // get the tree
	89	if (!fESDTree) {
	90	AliError("no ESD tree found");
	91	fESDFile->Close();
	92	fESDFile = 0x0;
	93	return;
	94	}
	95	fESDEvent->ReadFromTree(fESDTree); // link fESDEvent to the tree
068e4c36	96	}
068e4c36	97
7c29acf1	98	//_____________________________________________________________________________
	99	AliMUONRecoCheck::AliMUONRecoCheck(AliESDEvent esdEvent, AliMCEventHandler mcEventHandler)
	100	: TObject(),
	101	fMCEventHandler(0),
	102	fESDEvent(0),
	103	fESDTree (0x0),
	104	fESDFile (0x0),
	105	fCurrentEvent(0),
	106	fTrackRefStore(0x0),
	107	fRecoTrackRefStore(0x0),
	108	fRecoTrackStore(0x0),
	109	fESDEventOwner(kFALSE)
	110	{
	111	/// Normal ctor
	112
	113	// TrackRefs and Particules
	114	fMCEventHandler = mcEventHandler;
	115
	116	// ESD MUON Tracks
	117	fESDEvent = esdEvent;
	118
	119	}
	120
b8dc484b	121	//_____________________________________________________________________________
48217459	122	AliMUONRecoCheck::~AliMUONRecoCheck()
b8dc484b	123	{
48217459	124	/// Destructor
7c29acf1	125	if (fESDEventOwner) {
	126	delete fMCEventHandler;
	127	delete fESDEvent;
	128	if (fESDFile) fESDFile->Close();
	129	}
96ebe67e	130	ResetStores();
	131	}
	132
	133	//_____________________________________________________________________________
	134	void AliMUONRecoCheck::ResetStores()
	135	{
	136	/// Deletes all the store objects that have been created and resets the pointers to 0x0
	137	delete fTrackRefStore; fTrackRefStore = 0x0;
	138	delete fRecoTrackRefStore; fRecoTrackRefStore = 0x0;
	139	delete fRecoTrackStore; fRecoTrackStore = 0x0;
	140	}
	141
	142	//_____________________________________________________________________________
	143	Int_t AliMUONRecoCheck::NumberOfEvents() const
	144	{
	145	/// Return the number of events
7c29acf1	146	if (fESDEventOwner && fESDTree) return fESDTree->GetEntries();
96ebe67e	147	return 0;
b8dc484b	148	}
b8dc484b	149
48217459	150	//_____________________________________________________________________________
96ebe67e	151	AliMUONVTrackStore* AliMUONRecoCheck::ReconstructedTracks(Int_t event)
48217459	152	{
61fed964	153	/// Return a track store containing the reconstructed tracks (converted into
61fed964	154	/// MUONTrack objects) for a given event
7c29acf1	155
	156	if (!fESDEventOwner) {
	157	MakeReconstructedTracks();
	158	return fRecoTrackStore;
	159	}
	160
96ebe67e	161	if (event != fCurrentEvent) {
	162	ResetStores();
	163	fCurrentEvent = event;
	164	}
	165
	166	if (fRecoTrackStore != 0x0) return fRecoTrackStore;
61fed964	167	else {
	168	if (!fESDTree) return 0x0;
	169	if (fESDTree->GetEvent(event) <= 0) {
	170	AliError(Form("fails to read ESD object for event %d", event));
	171	return 0x0;
9bdda5f6	172	}
61fed964	173	MakeReconstructedTracks();
61fed964	174	return fRecoTrackStore;
9bdda5f6	175	}
9bdda5f6	176	}
	177
	178	//_____________________________________________________________________________
96ebe67e	179	AliMUONVTrackStore* AliMUONRecoCheck::TrackRefs(Int_t event)
b8dc484b	180	{
48217459	181	/// Return a track store containing the track references (converted into
48217459	182	/// MUONTrack objects) for a given event
7c29acf1	183
	184	if (!fESDEventOwner) {
	185	MakeTrackRefs();
	186	return fTrackRefStore;
	187	}
	188
96ebe67e	189	if (event != fCurrentEvent) {
	190	ResetStores();
	191	fCurrentEvent = event;
	192	}
	193
	194	if (fTrackRefStore != 0x0) return fTrackRefStore;
	195	else {
61fed964	196	if (!fMCEventHandler->GetEvent(event)) {
	197	AliError(Form("fails to read MC objects for event %d", event));
	198	return 0x0;
	199	}
96ebe67e	200	MakeTrackRefs();
	201	return fTrackRefStore;
	202	}
b8dc484b	203	}
	204
	205	//_____________________________________________________________________________
96ebe67e	206	AliMUONVTrackStore* AliMUONRecoCheck::ReconstructibleTracks(Int_t event)
b8dc484b	207	{
48217459	208	/// Return a track store containing the reconstructible tracks for a given event
7c29acf1	209
	210	if (!fESDEventOwner) {
	211	if (TrackRefs(event) == 0x0) return 0x0;
	212	MakeReconstructibleTracks();
	213	return fRecoTrackRefStore;
	214	}
	215
96ebe67e	216	if (event != fCurrentEvent) {
	217	ResetStores();
	218	fCurrentEvent = event;
	219	}
	220
	221	if (fRecoTrackRefStore != 0x0) return fRecoTrackRefStore;
	222	else {
	223	if (TrackRefs(event) == 0x0) return 0x0;
	224	MakeReconstructibleTracks();
	225	return fRecoTrackRefStore;
	226	}
48217459	227	}
b8dc484b	228
61fed964	229	//_____________________________________________________________________________
	230	void AliMUONRecoCheck::MakeReconstructedTracks()
	231	{
	232	/// Make reconstructed tracks
103e6575	233	if (!(fRecoTrackStore = AliMUONESDInterface::NewTrackStore())) return;
61fed964	234
b1fea02e	235	// loop over all reconstructed tracks and add them to the store (skip ghosts)
103e6575	236	Int_t nTracks = (Int_t) fESDEvent->GetNumberOfMuonTracks();
b1fea02e	237	for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
b1fea02e	238	AliESDMuonTrack* esdTrack = fESDEvent->GetMuonTrack(iTrack);
ad3c6eda	239	if (esdTrack->ContainTrackerData()) AliMUONESDInterface::Add(esdTrack, fRecoTrackStore);
b1fea02e	240	}
61fed964	241
	242	}
	243
48217459	244	//_____________________________________________________________________________
96ebe67e	245	void AliMUONRecoCheck::MakeTrackRefs()
48217459	246	{
48217459	247	/// Make reconstructible tracks
103e6575	248	AliMUONVTrackStore *tmpTrackRefStore = AliMUONESDInterface::NewTrackStore();
103e6575	249	if (!tmpTrackRefStore) return;
48217459	250
96ebe67e	251	Double_t x, y, z, pX, pY, pZ, bendingSlope, nonBendingSlope, inverseBendingMomentum;
	252	TParticle* particle;
	253	TClonesArray* trackRefs;
	254	Int_t nTrackRef = fMCEventHandler->MCEvent()->GetNumberOfTracks();
103e6575	255	AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
	256	if (!cStore) return;
	257	AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
b8dc484b	258
96ebe67e	259	// loop over simulated tracks
	260	for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; ++iTrackRef) {
	261	Int_t nHits = fMCEventHandler->GetParticleAndTR(iTrackRef, particle, trackRefs);
48217459	262
96ebe67e	263	// skip empty trackRefs
	264	if (nHits < 1) continue;
	265
61fed964	266	// get the particle charge for further calculation
61fed964	267	TParticlePDG* ppdg = particle->GetPDG();
46df197d	268	Int_t charge = ppdg != NULL ? (Int_t)(ppdg->Charge()/3.0) : 0;
96ebe67e	269
61fed964	270	AliMUONTrack track;
48217459	271
96ebe67e	272	// loop over simulated track hits
	273	for (Int_t iHit = 0; iHit < nHits; ++iHit) {
	274	AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(iHit));
068e4c36	275
61fed964	276	// skip trackRefs not in MUON
	277	if (trackReference->DetectorId() != AliTrackReference::kMUON) continue;
	278
96ebe67e	279	// Get track parameters of current hit
	280	x = trackReference->X();
	281	y = trackReference->Y();
	282	z = trackReference->Z();
	283	pX = trackReference->Px();
	284	pY = trackReference->Py();
	285	pZ = trackReference->Pz();
b8dc484b	286
96ebe67e	287	// check chamberId of current trackReference
c59f70b9	288	Int_t detElemId = trackReference->UserId();
c59f70b9	289	Int_t chamberId = detElemId / 100 - 1;
96ebe67e	290	if (chamberId < 0 \|\| chamberId >= AliMUONConstants::NTrackingCh()) continue;
48217459	291
96ebe67e	292	// set hit parameters
c59f70b9	293	hit->SetUniqueID(AliMUONVCluster::BuildUniqueID(chamberId, detElemId, iHit));
103e6575	294	hit->SetXYZ(x,y,z);
103e6575	295	hit->SetErrXY(0.,0.);
b8dc484b	296
96ebe67e	297	// compute track parameters at hit
7cdc832b	298	bendingSlope = 0;
	299	nonBendingSlope = 0;
	300	inverseBendingMomentum = 0;
96ebe67e	301	if (TMath::Abs(pZ) > 0) {
	302	bendingSlope = pY/pZ;
	303	nonBendingSlope = pX/pZ;
b8dc484b	304	}
96ebe67e	305	Double_t pYZ = TMath::Sqrt(pYpY+pZpZ);
	306	if (pYZ >0) inverseBendingMomentum = 1/pYZ;
	307	inverseBendingMomentum *= charge;
	308
	309	// set track parameters at hit
	310	AliMUONTrackParam trackParam;
	311	trackParam.SetNonBendingCoor(x);
	312	trackParam.SetBendingCoor(y);
	313	trackParam.SetZ(z);
	314	trackParam.SetBendingSlope(bendingSlope);
	315	trackParam.SetNonBendingSlope(nonBendingSlope);
	316	trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
	317
	318	// add track parameters at current hit to the track
103e6575	319	track.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
96ebe67e	320	}
96ebe67e	321
61fed964	322	// if none of the track hits was in MUON, goto the next track
	323	if (track.GetNClusters() < 1) continue;
	324
	325	// get track parameters at particle's vertex
	326	x = particle->Vx();
	327	y = particle->Vy();
	328	z = particle->Vz();
	329	pX = particle->Px();
	330	pY = particle->Py();
	331	pZ = particle->Pz();
	332
	333	// compute rest of track parameters at particle's vertex
	334	bendingSlope = 0;
	335	nonBendingSlope = 0;
	336	inverseBendingMomentum = 0;
	337	if (TMath::Abs(pZ) > 0) {
	338	bendingSlope = pY/pZ;
	339	nonBendingSlope = pX/pZ;
	340	}
	341	Double_t pYZ = TMath::Sqrt(pYpY+pZpZ);
	342	if (pYZ >0) inverseBendingMomentum = 1/pYZ;
	343	inverseBendingMomentum *= charge;
	344
	345	// set track parameters at particle's vertex
	346	AliMUONTrackParam trackParamAtVertex;
	347	trackParamAtVertex.SetNonBendingCoor(x);
	348	trackParamAtVertex.SetBendingCoor(y);
	349	trackParamAtVertex.SetZ(z);
	350	trackParamAtVertex.SetBendingSlope(bendingSlope);
	351	trackParamAtVertex.SetNonBendingSlope(nonBendingSlope);
	352	trackParamAtVertex.SetInverseBendingMomentum(inverseBendingMomentum);
	353
	354	// add track parameters at vertex
	355	track.SetTrackParamAtVertex(&trackParamAtVertex);
	356
103e6575	357	// store the track
2e2d0c44	358	track.SetUniqueID(iTrackRef);
96ebe67e	359	tmpTrackRefStore->Add(track);
b8dc484b	360	}
b8dc484b	361
96ebe67e	362	CleanMuonTrackRef(tmpTrackRefStore);
48217459	363
103e6575	364	delete hit;
103e6575	365	delete cStore;
96ebe67e	366	delete tmpTrackRefStore;
b8dc484b	367	}
	368
	369	//_____________________________________________________________________________
96ebe67e	370	void AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore *tmpTrackRefStore)
b8dc484b	371	{
48217459	372	/// Re-calculate hits parameters because two AliTrackReferences are recorded for
	373	/// each chamber (one when particle is entering + one when particle is leaving
	374	/// the sensitive volume)
103e6575	375	if (!(fTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
48217459	376
96ebe67e	377	Double_t maxGasGap = 1.; // cm
	378	Double_t x, y, z, pX, pY, pZ, x1, y1, z1, pX1, pY1, pZ1, z2;
	379	Double_t bendingSlope,nonBendingSlope,inverseBendingMomentum;
103e6575	380	AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
	381	if (!cStore) return;
	382	AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
48217459	383
96ebe67e	384	// create iterator
96ebe67e	385	TIter next(tmpTrackRefStore->CreateIterator());
48217459	386
96ebe67e	387	// loop over tmpTrackRef
48217459	388	AliMUONTrack* track;
96ebe67e	389	while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
	390
	391	AliMUONTrack newTrack;
	392
	393	// loop over tmpTrackRef's hits
	394	Int_t iHit1 = 0;
	395	Int_t nTrackHits = track->GetNClusters();
	396	while (iHit1 < nTrackHits) {
	397	AliMUONTrackParam trackParam1 = (AliMUONTrackParam) track->GetTrackParamAtCluster()->UncheckedAt(iHit1);
	398
	399	// get track parameters at hit1
	400	x1 = trackParam1->GetNonBendingCoor();
	401	y1 = trackParam1->GetBendingCoor();
	402	z1 = trackParam1->GetZ();
	403	pX1 = trackParam1->Px();
	404	pY1 = trackParam1->Py();
	405	pZ1 = trackParam1->Pz();
	406
	407	// prepare new track parameters
	408	x = x1;
	409	y = y1;
	410	z = z1;
	411	pX = pX1;
	412	pY = pY1;
	413	pZ = pZ1;
	414
	415	// loop over next tmpTrackRef's hits
	416	Int_t nCombinedHits = 1;
	417	for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
	418	AliMUONTrackParam trackParam2 = (AliMUONTrackParam) track->GetTrackParamAtCluster()->UncheckedAt(iHit2);
48217459	419
96ebe67e	420	// get z position of hit2
	421	z2 = trackParam2->GetZ();
	422
	423	// complete new track parameters if hit2 is on the same detection element
	424	if ( TMath::Abs(z2-z1) < maxGasGap ) {
	425	x += trackParam2->GetNonBendingCoor();
	426	y += trackParam2->GetBendingCoor();
	427	z += z2;
	428	pX += trackParam2->Px();
	429	pY += trackParam2->Py();
	430	pZ += trackParam2->Pz();
	431	nCombinedHits++;
48217459	432	iHit1 = iHit2;
	433	}
	434
96ebe67e	435	}
	436
	437	// finalize new track parameters
	438	x /= (Double_t)nCombinedHits;
	439	y /= (Double_t)nCombinedHits;
	440	z /= (Double_t)nCombinedHits;
	441	pX /= (Double_t)nCombinedHits;
	442	pY /= (Double_t)nCombinedHits;
	443	pZ /= (Double_t)nCombinedHits;
	444	bendingSlope = 0;
	445	nonBendingSlope = 0;
	446	inverseBendingMomentum = 0;
	447	if (TMath::Abs(pZ) > 0) {
	448	bendingSlope = pY/pZ;
	449	nonBendingSlope = pX/pZ;
	450	}
	451	Double_t pYZ = TMath::Sqrt(pYpY+pZpZ);
	452	if (pYZ >0) inverseBendingMomentum = 1/pYZ;
	453	inverseBendingMomentum *= trackParam1->GetCharge();
b8dc484b	454
96ebe67e	455	// set hit parameters
c59f70b9	456	hit->SetUniqueID(trackParam1->GetClusterPtr()->GetUniqueID());
103e6575	457	hit->SetXYZ(x,y,z);
103e6575	458	hit->SetErrXY(0.,0.);
48217459	459
96ebe67e	460	// set new track parameters at new hit
	461	AliMUONTrackParam trackParam;
	462	trackParam.SetNonBendingCoor(x);
	463	trackParam.SetBendingCoor(y);
	464	trackParam.SetZ(z);
	465	trackParam.SetBendingSlope(bendingSlope);
	466	trackParam.SetNonBendingSlope(nonBendingSlope);
	467	trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
	468
	469	// add track parameters at current hit to the track
103e6575	470	newTrack.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
b8dc484b	471
b8dc484b	472	iHit1++;
96ebe67e	473	}
48217459	474
2e2d0c44	475	newTrack.SetUniqueID(track->GetUniqueID());
96ebe67e	476	newTrack.SetTrackParamAtVertex(track->GetTrackParamAtVertex());
96ebe67e	477	fTrackRefStore->Add(newTrack);
b8dc484b	478
96ebe67e	479	}
b8dc484b	480
103e6575	481	delete hit;
103e6575	482	delete cStore;
b8dc484b	483	}
	484
	485	//_____________________________________________________________________________
96ebe67e	486	void AliMUONRecoCheck::MakeReconstructibleTracks()
b8dc484b	487	{
61fed964	488	/// Isolate the reconstructible tracks
103e6575	489	if (!(fRecoTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
48217459	490
96ebe67e	491	// create iterator on trackRef
96ebe67e	492	TIter next(fTrackRefStore->CreateIterator());
b8dc484b	493
96ebe67e	494	// loop over trackRef
48217459	495	AliMUONTrack* track;
96ebe67e	496	while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
	497
	498	Bool_t* chamberInTrack = new Bool_t(AliMUONConstants::NTrackingCh());
	499	for (Int_t iCh = 0; iCh < AliMUONConstants::NTrackingCh(); iCh++) chamberInTrack[iCh] = kFALSE;
48217459	500
96ebe67e	501	// loop over trackRef's hits to get hit chambers
	502	Int_t nTrackHits = track->GetNClusters();
	503	for (Int_t iHit = 0; iHit < nTrackHits; iHit++) {
	504	AliMUONVCluster* hit = ((AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit))->GetClusterPtr();
	505	chamberInTrack[hit->GetChamberId()] = kTRUE;
cc056de1	506	}
96ebe67e	507
cc056de1	508	// track is reconstructible if the particle is depositing a hit
cc056de1	509	// in the following chamber combinations:
96ebe67e	510	Bool_t trackOK = kTRUE;
	511	if (!chamberInTrack[0] && !chamberInTrack[1]) trackOK = kFALSE;
	512	if (!chamberInTrack[2] && !chamberInTrack[3]) trackOK = kFALSE;
	513	if (!chamberInTrack[4] && !chamberInTrack[5]) trackOK = kFALSE;
	514	Int_t nHitsInLastStations = 0;
	515	for (Int_t iCh = 6; iCh < AliMUONConstants::NTrackingCh(); iCh++)
	516	if (chamberInTrack[iCh]) nHitsInLastStations++;
	517	if(nHitsInLastStations < 3) trackOK = kFALSE;
48217459	518
96ebe67e	519	// Add reconstructible tracks to fRecoTrackRefStore
96ebe67e	520	if (trackOK) fRecoTrackRefStore->Add(*track);
48217459	521
96ebe67e	522	delete [] chamberInTrack;
b8dc484b	523	}
48217459	524
b8dc484b	525	}
b8dc484b	526