[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 "AliMUON.h"
#include "AliMUONRecoCheck.h"
#include "AliMUONTrack.h"
#include "AliMUONRecData.h"
#include "AliMUONConstants.h"

#include "AliLoader.h" 
#include "AliRunLoader.h" 
#include "AliHeader.h"
#include "AliStack.h"
#include "AliTrackReference.h"
#include "AliLog.h" 

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

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

//_____________________________________________________________________________
  AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader, Char_t *chLoaderSim)
  : TObject(),
  fRunLoader(0x0),
  fMUONData(0x0),
  fRunLoaderSim(0x0),
  fMUONDataSim(0x0),
  fMuonTrackRef(0x0),
  fTrackReco(0x0),
  fReconstructibleTracks(0),
  fRecoTracks(0),
  fIsLoadConstructor(kTRUE)
{
/// Constructor using "galice.root",
/// takes care of loading/unloading Kinematics, TrackRefs and Tracks

  fMuonTrackRef = new TClonesArray("AliMUONTrack", 10);

  // run loader
  fRunLoader = AliRunLoader::Open(chLoader,"MUONFolder","READ");
  if (!fRunLoader) {
    AliError(Form("no run loader found " ));
    return;
  }

 // initialize loader 	 
  AliLoader *loader = fRunLoader->GetLoader("MUONLoader");

  // container
  fMUONData  = new AliMUONRecData(loader,"MUON","MUON");
  if (!fMUONData) {
    AliError(Form("no MUONData found " ));
    return;
  }

  // run loader
  fRunLoaderSim = AliRunLoader::Open(chLoaderSim,"MUONFolderSim","READ");
  if (!fRunLoaderSim) {
    AliError(Form("no run sim loader found " ));
    return;
  }

 // initialize loader 	 
  AliLoader *loaderSim = fRunLoaderSim->GetLoader("MUONLoader");

  // container
  fMUONDataSim  = new AliMUONSimData(loaderSim,"MUON","MUON");
  if (!fMUONDataSim) {
    AliError(Form("no MUONDataSim found " ));
    return;
  }

  fRunLoaderSim->LoadKinematics("READ"); 	 
  fRunLoaderSim->LoadTrackRefs("READ"); 	 
  loader->LoadTracks("READ");

}

//_____________________________________________________________________________
  AliMUONRecoCheck::AliMUONRecoCheck(AliRunLoader *runloader, AliMUONRecData *muondata,
                                     AliRunLoader *runloaderSim, AliMUONSimData *muondataSim)
  : TObject(),
  fRunLoader(0x0),
  fMUONData(0x0),
  fRunLoaderSim(0x0),
  fMUONDataSim(0x0),
  fMuonTrackRef(0x0),
  fTrackReco(0x0),
  fReconstructibleTracks(0),
  fRecoTracks(0),
  fIsLoadConstructor(kFALSE)
{
/// Constructor from AliRunLoader and AliMUONData,
/// does not load/unload Kinematics, TrackRefs and Tracks internally
/// it should be done in the execution program (e.g. see MUONRecoCheck.C)

  fMuonTrackRef = new TClonesArray("AliMUONTrack", 10);

  // run loader
  fRunLoader = runloader;
  if (!fRunLoader) {
    AliError(Form("no run loader found " ));
    return;
  }

  // container
  fMUONData  = muondata;
  if (!fMUONData) {
    AliError(Form("no MUONData found " ));
    return;
  }

  // run loader
  fRunLoaderSim = runloaderSim;
  if (!fRunLoaderSim) {
    AliError(Form("no run sim loader found " ));
    return;
  }

  // container
  fMUONDataSim  = muondataSim;
  if (!fMUONDataSim) {
    AliError(Form("no MUONDataSim found " ));
    return;
  }

}


//_____________________________________________________________________________
AliMUONRecoCheck::~AliMUONRecoCheck()
{
/// Destructor
  delete fMuonTrackRef;
  if(fIsLoadConstructor){
    fRunLoaderSim->UnloadKinematics();
    fRunLoaderSim->UnloadTrackRefs();
    fRunLoader->UnloadTracks();
    delete fMUONData;
    delete fRunLoader;
    delete fMUONDataSim;
    delete fRunLoaderSim;
  }
}

//_____________________________________________________________________________
void AliMUONRecoCheck::MakeTrackRef()
{
/// Make reconstructible tracks

  AliTrackReference *trackReference;
  AliMUONTrack *muonTrack;
  AliMUONTrackParam *trackParam;
  AliMUONHitForRec *hitForRec;
  Float_t x, y, z, pX, pY, pZ, pYZ;
  Int_t track, trackSave;
  TParticle *particle;   
  Float_t bendingSlope = 0;
  Float_t  nonBendingSlope = 0;
  Float_t inverseBendingMomentum = 0;
  
  TTree* treeTR = fRunLoaderSim->TreeTR();
  if (treeTR == NULL) return;

  TBranch* branch = treeTR->GetBranch("MUON");
  if (branch == NULL) return;

  TClonesArray* trackRefs = 0;
  branch->SetAddress(&trackRefs);
  branch->SetAutoDelete(kTRUE);  

  Int_t nTrackRef = (Int_t)branch->GetEntries();
 
  track = trackSave = -999;
  Bool_t isNewTrack;
  Int_t iHitMin, iChamber, detElemId;

  trackParam = new AliMUONTrackParam();
  hitForRec = new AliMUONHitForRec();

  Int_t charge;
  TParticlePDG *ppdg;

  Int_t max = fRunLoaderSim->GetHeader()->Stack()->GetNtrack();
  for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; iTrackRef++) {

    branch->GetEntry(iTrackRef);
    
    iHitMin = 0;
    isNewTrack = kTRUE;
    
    if (!trackRefs->GetEntries()) continue; 

    while (isNewTrack) {

      muonTrack = new AliMUONTrack();
      
      for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); iHit++) {
      
	trackReference = (AliTrackReference*)trackRefs->At(iHit);
	x = trackReference->X();
	y = trackReference->Y();
	z = trackReference->Z();
	pX = trackReference->Px();
	pY = trackReference->Py();
	pZ = trackReference->Pz();

	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;
	}

	// track parameters at hit
	trackParam->SetBendingCoor(y);
	trackParam->SetNonBendingCoor(x);
	trackParam->SetZ(z);

	if (TMath::Abs(pZ) > 0) {
	  bendingSlope = pY/pZ;
	  nonBendingSlope = pX/pZ;
	}
	pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
	if (pYZ >0) inverseBendingMomentum = 1/pYZ; 

	trackParam->SetBendingSlope(bendingSlope);
	trackParam->SetNonBendingSlope(nonBendingSlope);
	trackParam->SetInverseBendingMomentum(inverseBendingMomentum);

	hitForRec->SetBendingCoor(y);
	hitForRec->SetNonBendingCoor(x);
	hitForRec->SetZ(z);
	hitForRec->SetBendingReso2(0.0); 
	hitForRec->SetNonBendingReso2(0.0);
	detElemId = hitForRec->GetDetElemId();
	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 
      particle = fRunLoaderSim->GetHeader()->Stack()->Particle(muonTrack->GetTrackID());

      if (particle) {

	x = particle->Vx();
	y = particle->Vy();
	z = particle->Vz();
	pX = particle->Px();
	pY = particle->Py();
	pZ = particle->Pz();

	trackParam->SetBendingCoor(y);
	trackParam->SetNonBendingCoor(x);
	trackParam->SetZ(z);
	if (TMath::Abs(pZ) > 0) {
	  bendingSlope = pY/pZ;
	  nonBendingSlope = pX/pZ;
	}

	pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
	if (pYZ >0) inverseBendingMomentum = 1/pYZ;       

	ppdg = particle->GetPDG(1);
	charge = (Int_t)(ppdg->Charge()/3.0);
	inverseBendingMomentum *= charge;

	trackParam->SetBendingSlope(bendingSlope);
	trackParam->SetNonBendingSlope(nonBendingSlope);
	trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
      
	muonTrack->SetTrackParamAtVertex(trackParam);
      }

      AddMuonTrackReference(muonTrack);
      delete muonTrack;
      muonTrack = NULL;
      
    } // end while isNewTrack
  }
  delete trackRefs;
  CleanMuonTrackRef();
  
  ReconstructibleTracks();

  delete trackParam;
  delete hitForRec;

}

//____________________________________________________________________________
TClonesArray* AliMUONRecoCheck::GetTrackReco()
{
/// Return TClonesArray of reconstructed tracks

  fMUONData->ResetRecTracks();
  fMUONData->SetTreeAddress("RT");
  fTrackReco = fMUONData->RecTracks(); 
  fMUONData->GetRecTracks();
  fRecoTracks = fTrackReco->GetEntriesFast();
  return fTrackReco;
}
//_____________________________________________________________________________
void AliMUONRecoCheck::PrintEvent() const
{
/// Debug facility

  AliMUONTrack *track;
  AliMUONHitForRec *hitForRec;
  TClonesArray *  hitForRecAtHit = 0;
  Float_t xRec,yRec,zRec;

  Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
  
  printf(" ******************************* \n");
  printf(" nb of tracks %d \n",nTrackRef);

  for (Int_t index = 0; index < nTrackRef; index++) {
    track = (AliMUONTrack*)fMuonTrackRef->At(index);
    hitForRecAtHit = track->GetHitForRecAtHit();
    Int_t nTrackHits = hitForRecAtHit->GetEntriesFast();
    printf(" track number %d \n",index);
    for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
      hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit); 
      xRec  = hitForRec->GetNonBendingCoor();
      yRec  = hitForRec->GetBendingCoor();
      zRec  = hitForRec->GetZ();
      printf("  x,y,z: %f , %f , %f \n",xRec,yRec,zRec);
    }
  }
}

//_____________________________________________________________________________
void AliMUONRecoCheck::ResetTracks() const
{
/// Reset tracks

  if (fMuonTrackRef) fMuonTrackRef->Delete();
}
//_____________________________________________________________________________
void AliMUONRecoCheck::CleanMuonTrackRef()
{
/// 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 
  AliMUONTrack *track, *trackNew;
  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;
  TClonesArray *newMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
  Int_t iHit1;
  Int_t iChamber = 0, detElemId = 0;
  Int_t nRec = 0;
  Int_t nTrackHits = 0;

  hitForRec = new AliMUONHitForRec();
  trackParam = new AliMUONTrackParam();

  Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
  for (Int_t index = 0; index < nTrackRef; index++) {
    track = (AliMUONTrack*)fMuonTrackRef->At(index);
    hitForRecAtHit = track->GetHitForRecAtHit();
    trackParamAtHit = track->GetTrackParamAtHit();
    trackParamAtVertex = track->GetTrackParamAtVertex();
    nTrackHits = hitForRecAtHit->GetEntriesFast();
    trackNew = new AliMUONTrack();
    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);
    {new ((*newMuonTrackRef)[newMuonTrackRef->GetEntriesFast()]) AliMUONTrack(*trackNew);}   
    delete trackNew;
    trackNew = NULL;
    
  } // end trackRef

  fMuonTrackRef->Delete();
  nTrackRef = newMuonTrackRef->GetEntriesFast();
  for (Int_t index = 0; index < nTrackRef; index++) {
    track = (AliMUONTrack*)newMuonTrackRef->At(index);
    AddMuonTrackReference(track);
  }
  
  delete hitForRec;
  delete trackParam;
  newMuonTrackRef->Delete();
  delete newMuonTrackRef;
  
}

//_____________________________________________________________________________
void AliMUONRecoCheck::ReconstructibleTracks()
{
/// Calculate the number of reconstructible tracks
  
  AliMUONTrack* track;
  TClonesArray* hitForRecAtHit = NULL;
  AliMUONHitForRec* hitForRec;
  Float_t zRec;
  Int_t nTrackRef, nTrackHits;
  Int_t isChamberInTrack[10];
  Int_t iChamber = 0;
  Bool_t isTrackOK = kTRUE;

  fReconstructibleTracks = 0;

  nTrackRef = fMuonTrackRef->GetEntriesFast();
  for (Int_t index = 0; index < nTrackRef; index++) {
    track = (AliMUONTrack*)fMuonTrackRef->At(index);
    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) fReconstructibleTracks++;
    if (!isTrackOK) fMuonTrackRef->Remove(track); // remove non reconstructible tracks
  }
  fMuonTrackRef->Compress();
}
Commit	Line	Data
b8dc484b	1	/**************************************************************************
	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	**************************************************************************/
	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 "AliMUON.h"
	25	#include "AliMUONRecoCheck.h"
	26	#include "AliMUONTrack.h"
6b092dfc	27	#include "AliMUONRecData.h"
cc056de1	28	#include "AliMUONConstants.h"
b8dc484b	29
b8dc484b	30	#include "AliLoader.h"
b8dc484b	31	#include "AliRunLoader.h"
b8dc484b	32	#include "AliHeader.h"
b8dc484b	33	#include "AliStack.h"
cc056de1	34	#include "AliTrackReference.h"
	35	#include "AliLog.h"
	36
	37	#include <TParticle.h>
5b1afea0	38	#include <TParticlePDG.h>
b8dc484b	39
5398f946	40	/// \cond CLASSIMP
b8dc484b	41	ClassImp(AliMUONRecoCheck)
5398f946	42	/// \endcond
b8dc484b	43
068e4c36	44	//_____________________________________________________________________________
6b092dfc	45	AliMUONRecoCheck::AliMUONRecoCheck(Char_t chLoader, Char_t chLoaderSim)
068e4c36	46	: TObject(),
	47	fRunLoader(0x0),
	48	fMUONData(0x0),
6b092dfc	49	fRunLoaderSim(0x0),
6b092dfc	50	fMUONDataSim(0x0),
068e4c36	51	fMuonTrackRef(0x0),
	52	fTrackReco(0x0),
	53	fReconstructibleTracks(0),
	54	fRecoTracks(0),
	55	fIsLoadConstructor(kTRUE)
	56	{
	57	/// Constructor using "galice.root",
	58	/// takes care of loading/unloading Kinematics, TrackRefs and Tracks
	59
	60	fMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
	61
	62	// run loader
6b092dfc	63	fRunLoader = AliRunLoader::Open(chLoader,"MUONFolder","READ");
068e4c36	64	if (!fRunLoader) {
	65	AliError(Form("no run loader found " ));
	66	return;
	67	}
	68
	69	// initialize loader
	70	AliLoader *loader = fRunLoader->GetLoader("MUONLoader");
	71
	72	// container
6b092dfc	73	fMUONData = new AliMUONRecData(loader,"MUON","MUON");
068e4c36	74	if (!fMUONData) {
	75	AliError(Form("no MUONData found " ));
	76	return;
	77	}
	78
6b092dfc	79	// run loader
	80	fRunLoaderSim = AliRunLoader::Open(chLoaderSim,"MUONFolderSim","READ");
	81	if (!fRunLoaderSim) {
	82	AliError(Form("no run sim loader found " ));
	83	return;
	84	}
	85
	86	// initialize loader
	87	AliLoader *loaderSim = fRunLoaderSim->GetLoader("MUONLoader");
	88
	89	// container
	90	fMUONDataSim = new AliMUONSimData(loaderSim,"MUON","MUON");
	91	if (!fMUONDataSim) {
	92	AliError(Form("no MUONDataSim found " ));
	93	return;
	94	}
	95
	96	fRunLoaderSim->LoadKinematics("READ");
cad8ba75	97	fRunLoaderSim->LoadTrackRefs("READ");
068e4c36	98	loader->LoadTracks("READ");
	99
	100	}
	101
b8dc484b	102	//_____________________________________________________________________________
6b092dfc	103	AliMUONRecoCheck::AliMUONRecoCheck(AliRunLoader runloader, AliMUONRecData muondata,
6b092dfc	104	AliRunLoader runloaderSim, AliMUONSimData muondataSim)
54d7ba50	105	: TObject(),
	106	fRunLoader(0x0),
	107	fMUONData(0x0),
6b092dfc	108	fRunLoaderSim(0x0),
6b092dfc	109	fMUONDataSim(0x0),
54d7ba50	110	fMuonTrackRef(0x0),
	111	fTrackReco(0x0),
	112	fReconstructibleTracks(0),
068e4c36	113	fRecoTracks(0),
068e4c36	114	fIsLoadConstructor(kFALSE)
b8dc484b	115	{
068e4c36	116	/// Constructor from AliRunLoader and AliMUONData,
	117	/// does not load/unload Kinematics, TrackRefs and Tracks internally
	118	/// it should be done in the execution program (e.g. see MUONRecoCheck.C)
5398f946	119
b8dc484b	120	fMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
b8dc484b	121
45a05f07	122	// run loader
45a05f07	123	fRunLoader = runloader;
b8dc484b	124	if (!fRunLoader) {
45a05f07	125	AliError(Form("no run loader found " ));
b8dc484b	126	return;
b8dc484b	127	}
b8dc484b	128
45a05f07	129	// container
	130	fMUONData = muondata;
	131	if (!fMUONData) {
	132	AliError(Form("no MUONData found " ));
	133	return;
	134	}
b8dc484b	135
6b092dfc	136	// run loader
	137	fRunLoaderSim = runloaderSim;
	138	if (!fRunLoaderSim) {
	139	AliError(Form("no run sim loader found " ));
	140	return;
	141	}
	142
	143	// container
	144	fMUONDataSim = muondataSim;
	145	if (!fMUONDataSim) {
	146	AliError(Form("no MUONDataSim found " ));
	147	return;
	148	}
	149
b8dc484b	150	}
b8dc484b	151
45a05f07	152
b8dc484b	153	//_____________________________________________________________________________
	154	AliMUONRecoCheck::~AliMUONRecoCheck()
	155	{
d19b6003	156	/// Destructor
b8dc484b	157	delete fMuonTrackRef;
068e4c36	158	if(fIsLoadConstructor){
6b092dfc	159	fRunLoaderSim->UnloadKinematics();
cad8ba75	160	fRunLoaderSim->UnloadTrackRefs();
068e4c36	161	fRunLoader->UnloadTracks();
068e4c36	162	delete fMUONData;
fac70e25	163	delete fRunLoader;
6b092dfc	164	delete fMUONDataSim;
6b092dfc	165	delete fRunLoaderSim;
068e4c36	166	}
b8dc484b	167	}
	168
	169	//_____________________________________________________________________________
	170	void AliMUONRecoCheck::MakeTrackRef()
	171	{
d19b6003	172	/// Make reconstructible tracks
b8dc484b	173
	174	AliTrackReference *trackReference;
	175	AliMUONTrack *muonTrack;
	176	AliMUONTrackParam *trackParam;
	177	AliMUONHitForRec *hitForRec;
	178	Float_t x, y, z, pX, pY, pZ, pYZ;
	179	Int_t track, trackSave;
	180	TParticle *particle;
	181	Float_t bendingSlope = 0;
	182	Float_t nonBendingSlope = 0;
	183	Float_t inverseBendingMomentum = 0;
	184
cad8ba75	185	TTree* treeTR = fRunLoaderSim->TreeTR();
b8dc484b	186	if (treeTR == NULL) return;
	187
	188	TBranch* branch = treeTR->GetBranch("MUON");
	189	if (branch == NULL) return;
	190
de876dd1	191	TClonesArray* trackRefs = 0;
b8dc484b	192	branch->SetAddress(&trackRefs);
de876dd1	193	branch->SetAutoDelete(kTRUE);
b8dc484b	194
de876dd1	195	Int_t nTrackRef = (Int_t)branch->GetEntries();
b8dc484b	196
	197	track = trackSave = -999;
	198	Bool_t isNewTrack;
30d3ea8a	199	Int_t iHitMin, iChamber, detElemId;
b8dc484b	200
	201	trackParam = new AliMUONTrackParam();
	202	hitForRec = new AliMUONHitForRec();
b8dc484b	203
5b1afea0	204	Int_t charge;
	205	TParticlePDG *ppdg;
	206
6b092dfc	207	Int_t max = fRunLoaderSim->GetHeader()->Stack()->GetNtrack();
b8dc484b	208	for (Int_t iTrackRef = 0; iTrackRef < nTrackRef; iTrackRef++) {
cc056de1	209
de876dd1	210	branch->GetEntry(iTrackRef);
b8dc484b	211
	212	iHitMin = 0;
	213	isNewTrack = kTRUE;
dba3c809	214
29fc2c86	215	if (!trackRefs->GetEntries()) continue;
b8dc484b	216
29fc2c86	217	while (isNewTrack) {
068e4c36	218
	219	muonTrack = new AliMUONTrack();
	220
b8dc484b	221	for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); iHit++) {
	222
	223	trackReference = (AliTrackReference*)trackRefs->At(iHit);
	224	x = trackReference->X();
	225	y = trackReference->Y();
	226	z = trackReference->Z();
	227	pX = trackReference->Px();
	228	pY = trackReference->Py();
	229	pZ = trackReference->Pz();
	230
	231	track = trackReference->GetTrack();
	232
cc056de1	233	if(track >= max){
	234	AliWarningStream()
	235	<< "Track ID " << track
	236	<< " larger than max number of particles " << max << endl;
	237	isNewTrack = kFALSE;
	238	break;
	239	}
b8dc484b	240	if (track != trackSave && iHit != 0) {
	241	iHitMin = iHit;
	242	trackSave = track;
	243	break;
	244	}
	245
	246	// track parameters at hit
	247	trackParam->SetBendingCoor(y);
	248	trackParam->SetNonBendingCoor(x);
	249	trackParam->SetZ(z);
cc056de1	250
b8dc484b	251	if (TMath::Abs(pZ) > 0) {
	252	bendingSlope = pY/pZ;
	253	nonBendingSlope = pX/pZ;
	254	}
	255	pYZ = TMath::Sqrt(pYpY+pZpZ);
	256	if (pYZ >0) inverseBendingMomentum = 1/pYZ;
	257
	258	trackParam->SetBendingSlope(bendingSlope);
	259	trackParam->SetNonBendingSlope(nonBendingSlope);
	260	trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
	261
	262	hitForRec->SetBendingCoor(y);
	263	hitForRec->SetNonBendingCoor(x);
	264	hitForRec->SetZ(z);
	265	hitForRec->SetBendingReso2(0.0);
30d3ea8a	266	hitForRec->SetNonBendingReso2(0.0);
	267	detElemId = hitForRec->GetDetElemId();
	268	if (detElemId) iChamber = detElemId / 100 - 1;
	269	else iChamber = AliMUONConstants::ChamberNumber(z);
b8dc484b	270	hitForRec->SetChamberNumber(iChamber);
b8dc484b	271
068e4c36	272	muonTrack->AddTrackParamAtHit(trackParam,0);
b8dc484b	273	muonTrack->AddHitForRecAtHit(hitForRec);
b8dc484b	274	muonTrack->SetTrackID(track);
cc056de1	275
b8dc484b	276	trackSave = track;
	277	if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
	278	}
cc056de1	279
b8dc484b	280	// track parameters at vertex
6b092dfc	281	particle = fRunLoaderSim->GetHeader()->Stack()->Particle(muonTrack->GetTrackID());
cc056de1	282
b8dc484b	283	if (particle) {
cc056de1	284
b8dc484b	285	x = particle->Vx();
	286	y = particle->Vy();
	287	z = particle->Vz();
	288	pX = particle->Px();
	289	pY = particle->Py();
	290	pZ = particle->Pz();
	291
	292	trackParam->SetBendingCoor(y);
c61c826b	293	trackParam->SetNonBendingCoor(x);
b8dc484b	294	trackParam->SetZ(z);
	295	if (TMath::Abs(pZ) > 0) {
	296	bendingSlope = pY/pZ;
	297	nonBendingSlope = pX/pZ;
	298	}
5b1afea0	299
b8dc484b	300	pYZ = TMath::Sqrt(pYpY+pZpZ);
b8dc484b	301	if (pYZ >0) inverseBendingMomentum = 1/pYZ;
5b1afea0	302
	303	ppdg = particle->GetPDG(1);
	304	charge = (Int_t)(ppdg->Charge()/3.0);
	305	inverseBendingMomentum *= charge;
	306
b8dc484b	307	trackParam->SetBendingSlope(bendingSlope);
	308	trackParam->SetNonBendingSlope(nonBendingSlope);
	309	trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
	310
	311	muonTrack->SetTrackParamAtVertex(trackParam);
	312	}
cc056de1	313
b8dc484b	314	AddMuonTrackReference(muonTrack);
068e4c36	315	delete muonTrack;
068e4c36	316	muonTrack = NULL;
b8dc484b	317
b8dc484b	318	} // end while isNewTrack
b8dc484b	319	}
2782bbd7	320	delete trackRefs;
b8dc484b	321	CleanMuonTrackRef();
	322
	323	ReconstructibleTracks();
	324
b8dc484b	325	delete trackParam;
b8dc484b	326	delete hitForRec;
b8dc484b	327
	328	}
	329
	330	//____________________________________________________________________________
	331	TClonesArray* AliMUONRecoCheck::GetTrackReco()
	332	{
d19b6003	333	/// Return TClonesArray of reconstructed tracks
b8dc484b	334
45a05f07	335	fMUONData->ResetRecTracks();
	336	fMUONData->SetTreeAddress("RT");
	337	fTrackReco = fMUONData->RecTracks();
	338	fMUONData->GetRecTracks();
b8dc484b	339	fRecoTracks = fTrackReco->GetEntriesFast();
	340	return fTrackReco;
	341	}
	342	//_____________________________________________________________________________
	343	void AliMUONRecoCheck::PrintEvent() const
	344	{
d19b6003	345	/// Debug facility
b8dc484b	346
	347	AliMUONTrack *track;
	348	AliMUONHitForRec *hitForRec;
	349	TClonesArray * hitForRecAtHit = 0;
	350	Float_t xRec,yRec,zRec;
	351
	352	Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
	353
	354	printf(" ******************************* \n");
	355	printf(" nb of tracks %d \n",nTrackRef);
	356
	357	for (Int_t index = 0; index < nTrackRef; index++) {
	358	track = (AliMUONTrack*)fMuonTrackRef->At(index);
	359	hitForRecAtHit = track->GetHitForRecAtHit();
	360	Int_t nTrackHits = hitForRecAtHit->GetEntriesFast();
	361	printf(" track number %d \n",index);
	362	for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
	363	hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
	364	xRec = hitForRec->GetNonBendingCoor();
	365	yRec = hitForRec->GetBendingCoor();
	366	zRec = hitForRec->GetZ();
	367	printf(" x,y,z: %f , %f , %f \n",xRec,yRec,zRec);
	368	}
	369	}
	370	}
	371
	372	//_____________________________________________________________________________
	373	void AliMUONRecoCheck::ResetTracks() const
	374	{
d19b6003	375	/// Reset tracks
d19b6003	376
2782bbd7	377	if (fMuonTrackRef) fMuonTrackRef->Delete();
b8dc484b	378	}
	379	//_____________________________________________________________________________
	380	void AliMUONRecoCheck::CleanMuonTrackRef()
	381	{
d19b6003	382	/// Re-calculate hits parameters because two AliTrackReferences are recorded for
	383	/// each chamber (one when particle is entering + one when particle is leaving
	384	/// the sensitive volume)
b8dc484b	385
	386	Float_t maxGasGap = 1.; // cm
	387	AliMUONTrack track, trackNew;
	388	AliMUONHitForRec hitForRec, hitForRec1, *hitForRec2;
	389	AliMUONTrackParam trackParam, trackParam1, trackParam2, trackParamAtVertex;
	390	TClonesArray * hitForRecAtHit = 0;
	391	TClonesArray * trackParamAtHit = 0;
	392	Float_t xRec,yRec,zRec;
	393	Float_t xRec1,yRec1,zRec1;
	394	Float_t xRec2,yRec2,zRec2;
	395	Float_t bendingSlope,nonBendingSlope,bendingMomentum;
	396	Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
	397	Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
	398	TClonesArray *newMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
	399	Int_t iHit1;
30d3ea8a	400	Int_t iChamber = 0, detElemId = 0;
b8dc484b	401	Int_t nRec = 0;
	402	Int_t nTrackHits = 0;
	403
	404	hitForRec = new AliMUONHitForRec();
	405	trackParam = new AliMUONTrackParam();
b8dc484b	406
b8dc484b	407	Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
b8dc484b	408	for (Int_t index = 0; index < nTrackRef; index++) {
	409	track = (AliMUONTrack*)fMuonTrackRef->At(index);
	410	hitForRecAtHit = track->GetHitForRecAtHit();
	411	trackParamAtHit = track->GetTrackParamAtHit();
	412	trackParamAtVertex = track->GetTrackParamAtVertex();
	413	nTrackHits = hitForRecAtHit->GetEntriesFast();
068e4c36	414	trackNew = new AliMUONTrack();
b8dc484b	415	iHit1 = 0;
	416	while (iHit1 < nTrackHits) {
	417	hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1);
	418	trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1);
	419	xRec1 = hitForRec1->GetNonBendingCoor();
	420	yRec1 = hitForRec1->GetBendingCoor();
	421	zRec1 = hitForRec1->GetZ();
	422	xRec = xRec1;
	423	yRec = yRec1;
	424	zRec = zRec1;
	425	bendingSlope1 = trackParam1->GetBendingSlope();
	426	nonBendingSlope1 = trackParam1->GetNonBendingSlope();
de876dd1	427	bendingMomentum1 = 0;
	428	if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
	429	bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
b8dc484b	430	bendingSlope = bendingSlope1;
	431	nonBendingSlope = nonBendingSlope1;
	432	bendingMomentum = bendingMomentum1;
	433	nRec = 1;
	434	for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
	435	hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2);
	436	trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2);
	437	xRec2 = hitForRec2->GetNonBendingCoor();
	438	yRec2 = hitForRec2->GetBendingCoor();
	439	zRec2 = hitForRec2->GetZ();
	440	bendingSlope2 = trackParam2->GetBendingSlope();
	441	nonBendingSlope2 = trackParam2->GetNonBendingSlope();
de876dd1	442	bendingMomentum2 = 0;
	443	if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
	444	bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
b8dc484b	445
b8dc484b	446	if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
cc056de1	447
b8dc484b	448	nRec++;
	449	xRec += xRec2;
	450	yRec += yRec2;
	451	zRec += zRec2;
	452	bendingSlope += bendingSlope2;
	453	nonBendingSlope += nonBendingSlope2;
	454	bendingMomentum += bendingMomentum2;
	455	iHit1 = iHit2;
	456	}
	457
	458	} // end iHit2
	459	xRec /= (Float_t)nRec;
	460	yRec /= (Float_t)nRec;
	461	zRec /= (Float_t)nRec;
	462	bendingSlope /= (Float_t)nRec;
	463	nonBendingSlope /= (Float_t)nRec;
	464	bendingMomentum /= (Float_t)nRec;
	465
	466	hitForRec->SetNonBendingCoor(xRec);
	467	hitForRec->SetBendingCoor(yRec);
	468	hitForRec->SetZ(zRec);
30d3ea8a	469	detElemId = hitForRec->GetDetElemId();
	470	if (detElemId) iChamber = detElemId / 100 - 1;
	471	else iChamber = AliMUONConstants::ChamberNumber(zRec);
b8dc484b	472	hitForRec->SetChamberNumber(iChamber);
	473	hitForRec->SetBendingReso2(0.0);
	474	hitForRec->SetNonBendingReso2(0.0);
	475	trackParam->SetNonBendingCoor(xRec);
	476	trackParam->SetBendingCoor(yRec);
	477	trackParam->SetZ(zRec);
	478	trackParam->SetNonBendingSlope(nonBendingSlope);
	479	trackParam->SetBendingSlope(bendingSlope);
de876dd1	480	if (TMath::Abs(bendingMomentum) > 0)
de876dd1	481	trackParam->SetInverseBendingMomentum(1./bendingMomentum);
b8dc484b	482
b8dc484b	483	trackNew->AddHitForRecAtHit(hitForRec);
068e4c36	484	trackNew->AddTrackParamAtHit(trackParam,0);
b8dc484b	485
	486	iHit1++;
	487	} // end iHit1
	488
	489	trackNew->SetTrackID(track->GetTrackID());
	490	trackNew->SetTrackParamAtVertex(trackParamAtVertex);
068e4c36	491	{new ((newMuonTrackRef)[newMuonTrackRef->GetEntriesFast()]) AliMUONTrack(trackNew);}
	492	delete trackNew;
	493	trackNew = NULL;
b8dc484b	494
	495	} // end trackRef
	496
2782bbd7	497	fMuonTrackRef->Delete();
b8dc484b	498	nTrackRef = newMuonTrackRef->GetEntriesFast();
	499	for (Int_t index = 0; index < nTrackRef; index++) {
	500	track = (AliMUONTrack*)newMuonTrackRef->At(index);
	501	AddMuonTrackReference(track);
	502	}
	503
b8dc484b	504	delete hitForRec;
b8dc484b	505	delete trackParam;
29fc2c86	506	newMuonTrackRef->Delete();
b8dc484b	507	delete newMuonTrackRef;
	508
	509	}
	510
	511	//_____________________________________________________________________________
	512	void AliMUONRecoCheck::ReconstructibleTracks()
	513	{
d19b6003	514	/// Calculate the number of reconstructible tracks
b8dc484b	515
	516	AliMUONTrack* track;
	517	TClonesArray* hitForRecAtHit = NULL;
	518	AliMUONHitForRec* hitForRec;
	519	Float_t zRec;
	520	Int_t nTrackRef, nTrackHits;
	521	Int_t isChamberInTrack[10];
	522	Int_t iChamber = 0;
	523	Bool_t isTrackOK = kTRUE;
	524
	525	fReconstructibleTracks = 0;
	526
	527	nTrackRef = fMuonTrackRef->GetEntriesFast();
	528	for (Int_t index = 0; index < nTrackRef; index++) {
	529	track = (AliMUONTrack*)fMuonTrackRef->At(index);
	530	hitForRecAtHit = track->GetHitForRecAtHit();
	531	nTrackHits = hitForRecAtHit->GetEntriesFast();
	532	for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
	533
	534	for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
	535	hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit);
	536	zRec = hitForRec->GetZ();
	537	iChamber = hitForRec->GetChamberNumber();
	538	if (iChamber < 0 \|\| iChamber > 10) continue;
	539	isChamberInTrack[iChamber] = 1;
cc056de1	540	}
	541	// track is reconstructible if the particle is depositing a hit
	542	// in the following chamber combinations:
	543
b8dc484b	544	isTrackOK = kTRUE;
cc056de1	545	if (!isChamberInTrack[0] && !isChamberInTrack[1]) isTrackOK = kFALSE;
	546	if (!isChamberInTrack[2] && !isChamberInTrack[3]) isTrackOK = kFALSE;
	547	if (!isChamberInTrack[4] && !isChamberInTrack[5]) isTrackOK = kFALSE;
	548	Int_t nHitsInLastStations=0;
	549	for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++)
	550	if (isChamberInTrack[ch]) nHitsInLastStations++;
	551	if(nHitsInLastStations < 3) isTrackOK = kFALSE;
	552
	553	if (isTrackOK) fReconstructibleTracks++;
b8dc484b	554	if (!isTrackOK) fMuonTrackRef->Remove(track); // remove non reconstructible tracks
	555	}
	556	fMuonTrackRef->Compress();
	557	}
	558