X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FAliMUONRecoCheck.cxx;h=5a83d2e54e687711da3c90390206a1f5a89eba7b;hb=12806ef7ba38a010261aea00f81c7758f8b391f3;hp=17b36c263f9b48e9faecf036d5a489ce224f0832;hpb=78649106ec7f71baeb897efa54274a4011675b87;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONRecoCheck.cxx b/MUON/AliMUONRecoCheck.cxx index 17b36c263f9..5a83d2e54e6 100644 --- a/MUON/AliMUONRecoCheck.cxx +++ b/MUON/AliMUONRecoCheck.cxx @@ -1,510 +1,526 @@ /************************************************************************** - * 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. * - **************************************************************************/ +* 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 +// 25 Jan 2008: +// Use the new ESDInterface to create MUON objects from ESD data +// - Philippe Pillot +// -#include "AliMUON.h" #include "AliMUONRecoCheck.h" #include "AliMUONTrack.h" -#include "AliMUONData.h" +#include "AliMUONVTrackStore.h" +#include "AliMUONVCluster.h" +#include "AliMUONVClusterStore.h" #include "AliMUONConstants.h" - -#include "AliLoader.h" -#include "AliRunLoader.h" -#include "AliHeader.h" +#include "AliMUONESDInterface.h" +#include "AliMUONTrackParam.h" +#include "AliMCEventHandler.h" +#include "AliMCEvent.h" #include "AliStack.h" #include "AliTrackReference.h" -#include "AliLog.h" +#include "AliLog.h" +#include "AliESDEvent.h" +#include "AliESDMuonTrack.h" +#include +#include #include +#include +#include /// \cond CLASSIMP ClassImp(AliMUONRecoCheck) /// \endcond //_____________________________________________________________________________ - AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader) - : TObject(), - fRunLoader(0x0), - fMUONData(0x0), - fMuonTrackRef(0x0), - fTrackReco(0x0), - fReconstructibleTracks(0), - fRecoTracks(0), - fIsLoadConstructor(kTRUE) +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) { -/// Constructor using "galice.root", -/// takes care of loading/unloading Kinematics, TrackRefs and Tracks - - fMuonTrackRef = new TClonesArray("AliMUONTrack", 10); - - // run loader - fRunLoader = AliRunLoader::Open(chLoader); - if (!fRunLoader) { - AliError(Form("no run loader found " )); + /// 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; } - - // initialize loader - AliLoader *loader = fRunLoader->GetLoader("MUONLoader"); - - // container - fMUONData = new AliMUONData(loader,"MUON","MUON"); - if (!fMUONData) { - AliError(Form("no MUONData found " )); + fESDTree = (TTree*) fESDFile->Get("esdTree"); // get the tree + if (!fESDTree) { + AliError("no ESD tree found"); + fESDFile->Close(); + fESDFile = 0x0; return; } - - fRunLoader->LoadKinematics("READ"); - fRunLoader->LoadTrackRefs("READ"); - loader->LoadTracks("READ"); - + fESDEvent->ReadFromTree(fESDTree); // link fESDEvent to the tree } //_____________________________________________________________________________ - AliMUONRecoCheck::AliMUONRecoCheck(AliRunLoader *runloader, AliMUONData *muondata) - : TObject(), - fRunLoader(0x0), - fMUONData(0x0), - fMuonTrackRef(0x0), - fTrackReco(0x0), - fReconstructibleTracks(0), - fRecoTracks(0), - fIsLoadConstructor(kFALSE) +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) { -/// 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; - } - + /// Normal ctor + + // TrackRefs and Particules + fMCEventHandler = mcEventHandler; + + // ESD MUON Tracks + fESDEvent = esdEvent; + } - //_____________________________________________________________________________ AliMUONRecoCheck::~AliMUONRecoCheck() { -/// Destructor - delete fMuonTrackRef; - if(fIsLoadConstructor){ - fRunLoader->UnloadKinematics(); - fRunLoader->UnloadTrackRefs(); - fRunLoader->UnloadTracks(); - delete fMUONData; + /// Destructor + if (fESDEventOwner) { + delete fMCEventHandler; + delete fESDEvent; + if (fESDFile) fESDFile->Close(); } - delete fRunLoader; + ResetStores(); } //_____________________________________________________________________________ -void AliMUONRecoCheck::MakeTrackRef() +void AliMUONRecoCheck::ResetStores() { -/// 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 = fRunLoader->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 max = fRunLoader->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 = fRunLoader->GetHeader()->Stack()->Particle(muonTrack->GetTrackID()); - - if (particle) { + /// 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; +} - x = particle->Vx(); - y = particle->Vy(); - z = particle->Vz(); - pX = particle->Px(); - pY = particle->Py(); - pZ = particle->Pz(); +//_____________________________________________________________________________ +Int_t AliMUONRecoCheck::NumberOfEvents() const +{ + /// Return the number of events + if (fESDEventOwner && fESDTree) return fESDTree->GetEntries(); + return 0; +} - 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); - - muonTrack->SetTrackParamAtVertex(trackParam); - } +//_____________________________________________________________________________ +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; + } - AddMuonTrackReference(muonTrack); - delete muonTrack; - muonTrack = NULL; - - } // end while isNewTrack + if (event != fCurrentEvent) { + ResetStores(); + fCurrentEvent = event; } - delete trackRefs; - CleanMuonTrackRef(); - ReconstructibleTracks(); - - delete trackParam; - delete hitForRec; - + 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; + } } -//____________________________________________________________________________ -TClonesArray* AliMUONRecoCheck::GetTrackReco() +//_____________________________________________________________________________ +AliMUONVTrackStore* AliMUONRecoCheck::TrackRefs(Int_t event) { -/// Return TClonesArray of reconstructed tracks + /// Return a track store containing the track references (converted into + /// MUONTrack objects) for a given event + + if (!fESDEventOwner) { + MakeTrackRefs(); + return fTrackRefStore; + } - fMUONData->ResetRecTracks(); - fMUONData->SetTreeAddress("RT"); - fTrackReco = fMUONData->RecTracks(); - fMUONData->GetRecTracks(); - fRecoTracks = fTrackReco->GetEntriesFast(); - return fTrackReco; + 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; + } } + //_____________________________________________________________________________ -void AliMUONRecoCheck::PrintEvent() const +AliMUONVTrackStore* AliMUONRecoCheck::ReconstructibleTracks(Int_t event) { -/// Debug facility + /// Return a track store containing the reconstructible tracks for a given event - AliMUONTrack *track; - AliMUONHitForRec *hitForRec; - TClonesArray * hitForRecAtHit = 0; - Float_t xRec,yRec,zRec; + if (!fESDEventOwner) { + if (TrackRefs(event) == 0x0) return 0x0; + MakeReconstructibleTracks(); + return fRecoTrackRefStore; + } - Int_t nTrackRef = fMuonTrackRef->GetEntriesFast(); + if (event != fCurrentEvent) { + ResetStores(); + fCurrentEvent = event; + } - 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); - } + if (fRecoTrackRefStore != 0x0) return fRecoTrackRefStore; + else { + if (TrackRefs(event) == 0x0) return 0x0; + MakeReconstructibleTracks(); + return fRecoTrackRefStore; } } //_____________________________________________________________________________ -void AliMUONRecoCheck::ResetTracks() const +void AliMUONRecoCheck::MakeReconstructedTracks() { -/// Reset tracks - - if (fMuonTrackRef) fMuonTrackRef->Delete(); + /// 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::CleanMuonTrackRef() +void AliMUONRecoCheck::MakeTrackRefs() { -/// 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(); + /// 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(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; +} - 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; +//_____________________________________________________________________________ +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(next()) ) ) { + + AliMUONTrack newTrack; + + // loop over tmpTrackRef's hits + Int_t iHit1 = 0; + Int_t nTrackHits = track->GetNClusters(); 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; + 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++) { - 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(); + AliMUONTrackParam *trackParam2 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit2); + + // get z position of hit2 + z2 = trackParam2->GetZ(); - 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; + // 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; + } + + } - 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); + // 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++; - } // end iHit1 - - trackNew->SetTrackID(track->GetTrackID()); - trackNew->SetTrackParamAtVertex(trackParamAtVertex); - {new ((*newMuonTrackRef)[newMuonTrackRef->GetEntriesFast()]) AliMUONTrack(*trackNew);} - delete trackNew; - trackNew = NULL; + } + + newTrack.SetUniqueID(track->GetUniqueID()); + newTrack.SetTrackParamAtVertex(track->GetTrackParamAtVertex()); + fTrackRefStore->Add(newTrack); - } // 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; - + delete hit; + delete cStore; } //_____________________________________________________________________________ -void AliMUONRecoCheck::ReconstructibleTracks() +void AliMUONRecoCheck::MakeReconstructibleTracks() { -/// Calculate the number of reconstructible tracks + /// Isolate the reconstructible tracks + if (!(fRecoTrackRefStore = AliMUONESDInterface::NewTrackStore())) return; + + // create iterator on trackRef + TIter next(fTrackRefStore->CreateIterator()); + // loop over trackRef 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; + while ( ( track = static_cast(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: - - 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 + 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; } - fMuonTrackRef->Compress(); + }