#if !defined(__CINT__) || defined(__MAKECINT__)
// ROOT includes
-#include <TTree.h>
#include <TFile.h>
#include <TH1.h>
#include <TCanvas.h>
#include <Riostream.h>
#include <TROOT.h>
#include <TClonesArray.h>
-#include <TGeoGlobalMagField.h>
+#include <TArrayI.h>
// STEER includes
#include "AliLog.h"
-#include "AliMagF.h"
#include "AliESDEvent.h"
#include "AliESDMuonTrack.h"
-#include "AliESDMuonCluster.h"
-#include "AliCDBPath.h"
-#include "AliCDBEntry.h"
#include "AliCDBManager.h"
// MUON includes
+#include "AliMUONCDB.h"
#include "AliMUONTrack.h"
#include "AliMUONVTrackStore.h"
#include "AliMUONTrackParam.h"
-#include "AliMUONTrackExtrap.h"
#include "AliMUONESDInterface.h"
#include "AliMUONRecoCheck.h"
#include "AliMUONVCluster.h"
/// Results are saved in the root file Fakes.root
/// Results are relevent provided that you use the same recoParams as for the reconstruction
-void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut);
-TTree* GetESDTree(TFile *esdFile);
-Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut);
-Bool_t IsRecontructible(AliMUONTrack &track, AliMUONRecoParam &recoParam);
-AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
- Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut);
-Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore, AliMUONRecoParam &recoParam,
- Bool_t useLabel, Double_t sigmaCut, TH1F &hFractionOfConnectedClusters);
+UInt_t requestedStationMask = 0;
+Bool_t request2ChInSameSt45 = kFALSE;
+Double_t sigmaCut = -1.;
+
+Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore,
+ Bool_t useLabel, TH1F &hFractionOfConnectedClusters);
//-----------------------------------------------------------------------
void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
- const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/")
+ const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/",
+ const TString ocdbPath = "local://$ALICE_ROOT/OCDB")
{
//Reset ROOT and connect tree file
TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.);
TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.);
- // prepare for analysis
- AliMUONRecoParam *recoParam = 0x0;
- Double_t sigmaCut = -1;
- Prepare(recoParam, sigmaCut);
+ // link to reconstructed and simulated tracks
+ AliMUONRecoCheck rc(esdFileName, SimDir);
- // link to reconstructed tracks
- TFile* esdFile = TFile::Open(esdFileName);
- TTree* esdTree = GetESDTree(esdFile);
- AliESDEvent* esd = new AliESDEvent();
- esd->ReadFromTree(esdTree);
+ // load necessary data from OCDB
+ AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
+ AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
+ AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
+ if (!recoParam) return;
- // link to simulated tracks
- AliMUONRecoCheck rc(esdFileName, SimDir);
+ // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
+ sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
+ // compute the mask of requested stations from recoParam
+ for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) requestedStationMask |= ( 1 << i );
+ // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
+ request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
// initialize global counters
+ Int_t nReconstructibleTracks = 0;
+ Int_t nReconstructedTracks = 0;
+ Int_t nEventsWithTrackReconstructedYet = 0;
Int_t nEventsWithFake = 0;
Int_t nEventsWithAdditionalFake = 0;
Int_t nTotMatchedTracks = 0;
Int_t nTotTracksReconstructedYet = 0;
Int_t nTotFakeTracks = 0;
+ Int_t nTotConnectedTracks = 0;
Int_t nTotAdditionalTracks = 0;
+ Bool_t trackReconstructedYet;
+ TArrayI eventsWithTrackReconstructedYet(10);
+ TArrayI eventsWithFake(10);
+ TArrayI eventsWithAdditionalFake(10);
// Loop over ESD events
FirstEvent = TMath::Max(0, FirstEvent);
- LastEvent = (LastEvent>=0) ? TMath::Min((Int_t)esdTree->GetEntries() - 1, LastEvent) : (Int_t)esdTree->GetEntries() - 1;
+ LastEvent = (LastEvent>=0) ? TMath::Min(rc.NumberOfEvents() - 1, LastEvent) : rc.NumberOfEvents() - 1;
for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
- // get the ESD of current event
- esdTree->GetEvent(iEvent);
- if (!esd) {
- Error("CheckESD", "no ESD object found for event %d", iEvent);
- return;
- }
-
- // convert TrackRef to MUON tracks
+ // get reconstructed and simulated tracks
+ AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(iEvent, kFALSE);
AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
+ if (!muonTrackStore || !trackRefStore) continue;
+
+ // count the number of reconstructible tracks
+ TIter next(trackRefStore->CreateIterator());
+ AliMUONTrack* trackRef;
+ while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
+ if (trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) nReconstructibleTracks++;
+ }
// loop over ESD tracks
Int_t nTrackerTracks = 0;
+ trackReconstructedYet = kFALSE;
AliMUONVTrackStore *fakeTrackStore = AliMUONESDInterface::NewTrackStore();
+ const AliESDEvent* esd = rc.GetESDEvent();
Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
- AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
+ AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
// skip ghosts
- if (!muonTrack->ContainTrackerData()) continue;
+ if (!esdTrack->ContainTrackerData()) continue;
nTrackerTracks++;
+ // find the corresponding MUON track
+ AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID());
+
// get track info
- Int_t nClusters = muonTrack->GetNClusters();
- Double_t normalizedChi2 = muonTrack->GetChi2() / (2. * muonTrack->GetNHit() - 5);
- Double_t p = muonTrack->P();
- Double_t pT = muonTrack->Pt();
- Double_t eta = muonTrack->Eta();
- Double_t phi = muonTrack->Phi();
+ Int_t nClusters = esdTrack->GetNClusters();
+ Double_t normalizedChi2 = esdTrack->GetChi2() / (2. * esdTrack->GetNHit() - 5);
+ Double_t p = esdTrack->P();
+ Double_t pT = esdTrack->Pt();
+ Double_t eta = esdTrack->Eta();
+ Double_t phi = esdTrack->Phi();
// fill global histograms
hNumberOfClusters->Fill(nClusters);
hPhi->Fill(phi);
// try to match the reconstructed track with a simulated one
- Float_t fractionOfMatchCluster = 0.;
- AliMUONTrack* matchedTrackRef = MatchWithTrackRef(*muonTrack, *trackRefStore, fractionOfMatchCluster, useLabel, sigmaCut);
+ Int_t nMatchClusters = 0;
+ AliMUONTrack* matchedTrackRef = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClusters, useLabel, sigmaCut);
// take actions according to matching result
if (matchedTrackRef) {
// global counter
nTotMatchedTracks++;
- if (!IsRecontructible(*matchedTrackRef,*recoParam)) nTotTracksReconstructedYet++;
+ if (!matchedTrackRef->IsValid(requestedStationMask, request2ChInSameSt45)) {
+ trackReconstructedYet = kTRUE;
+ nTotTracksReconstructedYet++;
+ }
// fill histograms
- hFractionOfMatchedClusters->Fill(fractionOfMatchCluster);
+ hFractionOfMatchedClusters->Fill(((Float_t) nMatchClusters) / ((Float_t) nClusters));
hNumberOfClustersMC->Fill(matchedTrackRef->GetNClusters());
hNumberOfClustersM->Fill(nClusters);
hChi2PerDofM->Fill(normalizedChi2);
hPhiF->Fill(phi);
// store fake tracks
- AliMUONESDInterface::Add(*muonTrack, *fakeTrackStore);
+ fakeTrackStore->Add(*muonTrack);
}
// fill histograms
hNumberOfTracks->Fill(nTrackerTracks);
+ nReconstructedTracks += nTrackerTracks;
+ if (trackReconstructedYet) eventsWithTrackReconstructedYet[nEventsWithTrackReconstructedYet++] = iEvent;
// count the number the additional fake tracks
if (fakeTrackStore->GetSize() > 0) {
// remove the most connected fake tracks
- Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore, *recoParam,
- useLabel, sigmaCut, *hFractionOfConnectedClusters);
+ Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore, useLabel, *hFractionOfConnectedClusters);
// remove the remaining free reconstructible tracks
Int_t nAdditionalTracks = fakeTrackStore->GetSize() - nFreeMissingTracks;
// fill histograms
+ eventsWithFake[nEventsWithFake] = iEvent;
nEventsWithFake++;
if (nAdditionalTracks > 0) {
+ eventsWithAdditionalFake[nEventsWithAdditionalFake] = iEvent;
nEventsWithAdditionalFake++;
nTotAdditionalTracks += nAdditionalTracks;
hNumberOfAdditionalTracks->Fill(nAdditionalTracks);
delete fakeTrackStore;
} // end of loop over events
+
+ // total number of connected tracks
+ nTotConnectedTracks = hFractionOfConnectedClusters->GetEntries();
// plot results
TCanvas cFakesSummary("cFakesSummary","cFakesSummary",900,600);
// print results
cout << endl;
+ cout << "- Number of reconstructible tracks: " << nReconstructibleTracks << endl;
+ cout << "- Number of reconstructed tracks: " << nReconstructedTracks << endl;
cout << "- Number of matched tracks: " << nTotMatchedTracks << endl;
- cout << " (including " << nTotTracksReconstructedYet << " tracks matched with a TrackRef that is not reconstructible)" << endl;
+ cout << " (including " << nTotTracksReconstructedYet << " track(s) matched with a TrackRef that is not reconstructible";
+ if (nTotTracksReconstructedYet > 0) {
+ for(Int_t i=0; i<nEventsWithTrackReconstructedYet; i++){
+ if (i==0) cout << " (eventID = " << eventsWithTrackReconstructedYet[i];
+ else cout << ", " << eventsWithTrackReconstructedYet[i];
+ }
+ cout << "))" << endl;
+ } else cout << ")" << endl;
cout << "- Number of fake tracks: " << nTotFakeTracks << endl;
- cout << " (including " << nTotAdditionalTracks << " additional tracks (compared to the number of expected ones))" << endl;
- cout << "- Number of events with fake track(s): " << nEventsWithFake << endl;
- cout << " (including " << nEventsWithAdditionalFake << " events with additional tracks)" << endl;
+ cout << " (including " << nTotConnectedTracks << " track(s) still connected to a reconstructible one)" << endl;
+ cout << " (including " << nTotAdditionalTracks << " additional track(s) (compared to the number of expected ones))" << endl;
+ cout << "- Number of events with fake track(s): " << nEventsWithFake;
+ if (nEventsWithFake > 0) {
+ for(Int_t i=0; i<nEventsWithFake; i++){
+ if (i==0) cout << " (eventID = " << eventsWithFake[i];
+ else cout << ", " << eventsWithFake[i];
+ }
+ cout << ")" << endl;
+ } else cout << endl;
+ cout << " (including " << nEventsWithAdditionalFake << " events with additional track(s)";
+ if (nEventsWithAdditionalFake > 0) {
+ for(Int_t i=0; i<nEventsWithAdditionalFake; i++){
+ if (i==0) cout << " (eventID = " << eventsWithAdditionalFake[i];
+ else cout << ", " << eventsWithAdditionalFake[i];
+ }
+ cout << "))" << endl;
+ } else cout << ")" << endl;
cout << endl;
cout << "REMINDER: results are relevent provided that you use the same recoParams as for the reconstruction" << endl;
cout << endl;
- // clear memory
- delete esd;
- esdFile->Close();
- delete recoParam;
-
-}
-
-//-----------------------------------------------------------------------
-void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut)
-{
- /// Set the magnetic field and return recoParam and sigmaCut to associate cluster and trackRef
-
- // prepare OCDB access
- AliCDBManager* man = AliCDBManager::Instance();
- man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
- man->SetRun(0);
-
- // set mag field
- // waiting for mag field in CDB
- if (!TGeoGlobalMagField::Instance()->GetField()) {
- printf("Loading field map...\n");
- AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG);
- TGeoGlobalMagField::Instance()->SetField(field);
- }
- // set the magnetic field for track extrapolations
- AliMUONTrackExtrap::SetField();
-
- // Load initial reconstruction parameters from OCDB
- AliCDBPath path("MUON","Calib","RecoParam");
- AliCDBEntry *entry=man->Get(path.GetPath());
- if(entry) {
- recoParam = dynamic_cast<AliMUONRecoParam*>(entry->GetObject());
- entry->SetOwner(0);
- AliCDBManager::Instance()->UnloadFromCache(path.GetPath());
- }
- if (!recoParam) {
- printf("Couldn't find RecoParam object in OCDB: create default one");
- recoParam = AliMUONRecoParam::GetLowFluxParam();
- }
-
- Info("MUONFakes", "\n recontruction parameters:");
- recoParam->Print("FULL");
- AliMUONESDInterface::ResetTracker(recoParam);
-
- // sigma cut to associate clusters with TrackRefs in case the label are not used
- sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
-
-}
-
-//-----------------------------------------------------------------------
-TTree* GetESDTree(TFile *esdFile)
-{
- /// Check that the file is properly open
- /// Return pointer to the ESD Tree
-
- if (!esdFile || !esdFile->IsOpen()) {
- Error("GetESDTree", "opening ESD file failed");
- exit(-1);
- }
-
- TTree* tree = (TTree*) esdFile->Get("esdTree");
- if (!tree) {
- Error("GetESDTree", "no ESD tree found");
- exit(-1);
- }
-
- return tree;
-
-}
-
-//-----------------------------------------------------------------------
-Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut)
-{
- /// Try to match 2 tracks
-
- Bool_t compTrack[10];
- Int_t nMatchClusters = track.CompatibleTrack(&trackRef, sigmaCut, compTrack);
- fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)track.GetNClusters());
-
- if ((compTrack[0] || compTrack[1] || compTrack[2] || compTrack[3]) && // at least 1 cluster matched in st 1 & 2
- (compTrack[6] || compTrack[7] || compTrack[8] || compTrack[9]) && // at least 1 cluster matched in st 4 & 5
- fractionOfMatchCluster > 0.5) return kTRUE; // more than 50% of clusters matched
- else return kFALSE;
-
-}
-
-//-----------------------------------------------------------------------
-Bool_t IsRecontructible(AliMUONTrack &track, AliMUONRecoParam &recoParam)
-{
- /// Check il the track is reconstructible
- Int_t nMinChHitInSt45 = (recoParam.MakeMoreTrackCandidates()) ? 2 : 3;
- Int_t currentCh, previousCh = -1, nChHitInSt45 = 0;
- Bool_t clusterInSt[5];
- for (Int_t iSt = 0; iSt < 5; iSt++) clusterInSt[iSt] = !recoParam.RequestStation(iSt);
-
- AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First());
- while (trackParam) {
-
- currentCh = trackParam->GetClusterPtr()->GetChamberId();
-
- clusterInSt[currentCh/2] = kTRUE;
-
- if (currentCh > 5 && currentCh != previousCh) {
- nChHitInSt45++;
- previousCh = currentCh;
- }
-
- trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
- }
-
- return (clusterInSt[0] && clusterInSt[1] && clusterInSt[2] &&
- clusterInSt[3] && clusterInSt[4] && nChHitInSt45 >= nMinChHitInSt45);
-
-}
-
-//-----------------------------------------------------------------------
-AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
- Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut)
-{
- /// Return if the trackRef matched with the reconstructed track and the fraction of matched clusters
-
- AliMUONTrack *matchedTrackRef = 0x0;
- fractionOfMatchCluster = 0.;
-
- if (useLabel) { // by using the MC label
-
- // get the corresponding simulated track if any
- Int_t label = muonTrack.GetLabel();
- matchedTrackRef = (AliMUONTrack*) trackRefStore.FindObject(label);
-
- // get the fraction of matched clusters
- if (matchedTrackRef) {
- Int_t nMatchClusters = 0;
- if (muonTrack.ClustersStored()) {
- AliESDMuonCluster* cluster = (AliESDMuonCluster*) muonTrack.GetClusters().First();
- while (cluster) {
- if (cluster->GetLabel() == label) nMatchClusters++;
- cluster = (AliESDMuonCluster*) muonTrack.GetClusters().After(cluster);
- }
- }
- fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)muonTrack.GetNClusters());
- }
-
- } else { // by comparing cluster/TrackRef positions
-
- // convert ESD track to MUON track
- AliMUONTrack track;
- AliMUONESDInterface::ESDToMUON(muonTrack,track);
-
- // look for the corresponding simulated track if any
- TIter next(trackRefStore.CreateIterator());
- AliMUONTrack* trackRef;
- while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
-
- // check compatibility
- Float_t f = 0.;
- if (TrackMatched(track, *trackRef, f, sigmaCut)) {
- matchedTrackRef = trackRef;
- fractionOfMatchCluster = f;
- break;
- }
-
- }
-
- }
-
- return matchedTrackRef;
-
}
//-----------------------------------------------------------------------
-Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore, AliMUONRecoParam &recoParam,
- Bool_t useLabel, Double_t sigmaCut, TH1F &hFractionOfConnectedClusters)
+Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore,
+ Bool_t useLabel, TH1F &hFractionOfConnectedClusters)
{
/// loop over reconstructible TrackRef not associated with reconstructed track:
/// for each of them, find and remove the most connected the fake track, if any,
while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
// skip not reconstructible trackRefs
- if (!IsRecontructible(*trackRef,recoParam)) continue;
+ if (!trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) continue;
Int_t label = trackRef->GetUniqueID();
// look for the most connected fake track
AliMUONTrack *connectedFake = 0x0;
- Float_t fractionOfConnectedClusters = 0.;
+ Double_t fractionOfConnectedClusters = 0.;
TIter next2(fakeTrackStore.CreateIterator());
AliMUONTrack* fakeTrack;
while ( ( fakeTrack = static_cast<AliMUONTrack*>(next2()) ) ) {
nConnectedClusters++;
} else { // by comparing cluster/TrackRef positions
Bool_t compTrack[10];
- nConnectedClusters = fakeTrack->CompatibleTrack(trackRef, sigmaCut, compTrack);
+ nConnectedClusters = fakeTrack->FindCompatibleClusters(*trackRef, sigmaCut, compTrack);
}
// skip non-connected fake tracks
if (nConnectedClusters == 0) continue;
// check if it is the most connected fake track
- Float_t f = ((Float_t)nConnectedClusters) / ((Float_t)fakeTrack->GetNClusters());
+ Double_t f = ((Double_t)nConnectedClusters) / ((Double_t)fakeTrack->GetNClusters());
if (f > fractionOfConnectedClusters) {
connectedFake = fakeTrack;
fractionOfConnectedClusters = f;
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff