1 #if !defined(__CINT__) || defined(__MAKECINT__)
8 #include <TClonesArray.h>
9 #include <TLorentzVector.h>
13 #include "AliESDEvent.h"
14 #include "AliESDMuonTrack.h"
15 #include "AliCDBManager.h"
18 #include "AliMUONCDB.h"
19 #include "AliMUONTrack.h"
20 #include "AliMUONVTrackStore.h"
21 #include "AliMUONTrackParam.h"
22 #include "AliMUONRecoCheck.h"
23 #include "AliMUONVCluster.h"
24 #include "AliMUONRecoParam.h"
28 /// \file DIMUONFakes.C
30 /// \author Ph. Pillot, Subatech, March. 2009
32 /// Macro to study the effects of fake tracks on the dimuon spectra
33 /// Results are saved in the root file DiFakes.root
34 /// Results are relevent provided that you use the same recoParams as for the reconstruction
36 Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut);
37 AliMUONTrack* MatchWithTrackRef(AliMUONTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
38 Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut);
40 //-----------------------------------------------------------------------
41 void DIMUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
42 const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/",
43 const TString ocdbPath = "local://$ALICE_ROOT/OCDB")
46 //Reset ROOT and connect tree file
49 // File for histograms and histogram booking
50 TFile *histoFile = new TFile("DiFakes.root", "RECREATE");
52 TH1F *hMass = new TH1F("hMass", "Dimuon mass distribution (GeV/c^{2})", 100, 0., 12.);
53 TH1F *hMassM = new TH1F("hMassM", "matched track mass distribution (GeV/c^{2})", 100, 0., 12.);
54 TH1F *hMassF = new TH1F("hMassF", "fake track mass distribution (GeV/c^{2})", 100, 0., 12.);
55 TH1F *hP = new TH1F("hP", "Dimuon P distribution (GeV/c)", 100, 0., 200.);
56 TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
57 TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
58 TH1F *hPt = new TH1F("hPt", "Dimuon Pt distribution (GeV/c)", 100, 0., 20.);
59 TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
60 TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
61 TH1F *hY = new TH1F("hY"," Dimuon rapidity distribution",100,-10,0);
62 TH1F *hYM = new TH1F("hYM"," matched track rapidity distribution",100,-10,0);
63 TH1F *hYF = new TH1F("hYF"," fake track rapidity distribution",100,-10,0);
64 TH1F *hEta = new TH1F("hEta"," Dimuon pseudo-rapidity distribution",100,-10,0);
65 TH1F *hEtaM = new TH1F("hEtaM"," matched track pseudo-rapidity distribution",100,-10,0);
66 TH1F *hEtaF = new TH1F("hEtaF"," fake track pseudo-rapidity distribution",100,-10,0);
67 TH1F *hPhi = new TH1F("hPhi"," Dimuon phi distribution",100,-1.,9.);
68 TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.);
69 TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.);
71 // link to reconstructed and simulated tracks
72 AliMUONRecoCheck rc(esdFileName, SimDir);
74 // load necessary data from OCDB
75 AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
76 AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
77 AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
78 if (!recoParam) return;
80 // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
81 Double_t sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
83 TLorentzVector vMu1, vMu2, vDiMu;
85 // Loop over ESD events
86 FirstEvent = TMath::Max(0, FirstEvent);
87 LastEvent = (LastEvent>=0) ? TMath::Min(rc.NumberOfEvents() - 1, LastEvent) : rc.NumberOfEvents() - 1;
88 for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
90 // get reconstructed and simulated tracks
91 AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(iEvent, kFALSE);
92 AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
93 if (!muonTrackStore || !trackRefStore) continue;
95 // loop over ESD tracks and flag them
96 const AliESDEvent* esd = rc.GetESDEvent();
97 Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
98 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
100 AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
103 if (!esdTrack->ContainTrackerData()) continue;
105 // find the corresponding MUON track
106 AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID());
108 // try to match the reconstructed track with a simulated one
109 Float_t fractionOfMatchCluster = 0.;
110 AliMUONTrack* matchedTrackRef = MatchWithTrackRef(*muonTrack, *trackRefStore, fractionOfMatchCluster, useLabel, sigmaCut);
112 // take actions according to matching result
113 if (matchedTrackRef) {
115 // flag matched tracks
116 esdTrack->SetLabel(matchedTrackRef->GetUniqueID());
118 // remove already matched trackRefs
119 trackRefStore->Remove(*matchedTrackRef);
124 esdTrack->SetLabel(-1);
130 // double loop over ESD tracks, build pairs and fill histograms according to their label
131 for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
132 AliESDMuonTrack* muonTrack1 = esd->GetMuonTrack(iTrack1);
135 if (!muonTrack1->ContainTrackerData()) continue;
138 Short_t charge1 = muonTrack1->Charge();
139 Int_t label1 = muonTrack1->GetLabel();
140 muonTrack1->LorentzP(vMu1);
142 for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
143 AliESDMuonTrack* muonTrack2 = esd->GetMuonTrack(iTrack2);
146 if (!muonTrack2->ContainTrackerData()) continue;
148 // keep only opposite sign pairs
149 Short_t charge2 = muonTrack2->Charge();
150 if (charge1*charge2 > 0) continue;
153 Int_t label2 = muonTrack2->GetLabel();
154 muonTrack2->LorentzP(vMu2);
156 // compute kinematics of the pair
158 Float_t mass = vDiMu.M();
159 Float_t p = vDiMu.P();
160 Float_t pt = vDiMu.Pt();
161 Float_t y = vDiMu.Rapidity();
162 Float_t eta = vDiMu.Eta();
163 Float_t phi = vDiMu.Phi();
164 if (phi < 0) phi += 2.*TMath::Pi();
166 // fill global histograms
174 // fill histograms according to labels
175 if (label1 >= 0 && label2 >= 0) {
199 } // end of loop over events
202 TCanvas cDiFakesSummary("cDiFakesSummary","cDiFakesSummary",900,600);
203 cDiFakesSummary.Divide(3,2);
204 cDiFakesSummary.cd(1);
205 cDiFakesSummary.GetPad(1)->SetLogy();
207 hMass->SetMinimum(0.5);
208 hMassM->Draw("same");
209 hMassM->SetLineColor(4);
210 hMassF->Draw("same");
211 hMassF->SetLineColor(2);
212 hMassF->SetFillColor(2);
213 hMassF->SetFillStyle(3017);
214 cDiFakesSummary.cd(2);
215 cDiFakesSummary.GetPad(3)->SetLogy();
219 hPM->SetLineColor(4);
221 hPF->SetLineColor(2);
222 hPF->SetFillColor(2);
223 hPF->SetFillStyle(3017);
224 cDiFakesSummary.cd(3);
225 cDiFakesSummary.GetPad(4)->SetLogy();
227 hPt->SetMinimum(0.5);
229 hPtM->SetLineColor(4);
231 hPtF->SetLineColor(2);
232 hPtF->SetFillColor(2);
233 hPtF->SetFillStyle(3017);
234 cDiFakesSummary.cd(4);
235 cDiFakesSummary.GetPad(2)->SetLogy();
239 hYM->SetLineColor(4);
241 hYF->SetLineColor(2);
242 hYF->SetFillColor(2);
243 hYF->SetFillStyle(3017);
244 cDiFakesSummary.cd(5);
245 cDiFakesSummary.GetPad(5)->SetLogy();
247 hEta->SetMinimum(0.5);
249 hEtaM->SetLineColor(4);
251 hEtaF->SetLineColor(2);
252 hEtaF->SetFillColor(2);
253 hEtaF->SetFillStyle(3017);
254 cDiFakesSummary.cd(6);
255 cDiFakesSummary.GetPad(6)->SetLogy();
257 hPhi->SetMinimum(0.5);
259 hPhiM->SetLineColor(4);
261 hPhiF->SetLineColor(2);
262 hPhiF->SetFillColor(2);
263 hPhiF->SetFillStyle(3017);
268 cDiFakesSummary.Write();
273 //-----------------------------------------------------------------------
274 Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut)
276 /// Try to match 2 tracks
278 Bool_t compTrack[10];
279 Int_t nMatchClusters = track.CompatibleTrack(&trackRef, sigmaCut, compTrack);
280 fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)track.GetNClusters());
282 if ((compTrack[0] || compTrack[1] || compTrack[2] || compTrack[3]) && // at least 1 cluster matched in st 1 & 2
283 (compTrack[6] || compTrack[7] || compTrack[8] || compTrack[9]) && // at least 1 cluster matched in st 4 & 5
284 fractionOfMatchCluster > 0.5) return kTRUE; // more than 50% of clusters matched
289 //-----------------------------------------------------------------------
290 AliMUONTrack* MatchWithTrackRef(AliMUONTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
291 Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut)
293 /// Return if the trackRef matched with the reconstructed track and the fraction of matched clusters
295 AliMUONTrack *matchedTrackRef = 0x0;
296 fractionOfMatchCluster = 0.;
298 if (useLabel) { // by using the MC label
300 // get the corresponding simulated track if any
301 Int_t label = muonTrack.GetMCLabel();
302 matchedTrackRef = (AliMUONTrack*) trackRefStore.FindObject(label);
304 // get the fraction of matched clusters
305 if (matchedTrackRef) {
306 Int_t nMatchClusters = 0;
307 Int_t nClusters = muonTrack.GetNClusters();
308 for (Int_t iCl = 0; iCl < nClusters; iCl++)
309 if (((AliMUONTrackParam*) muonTrack.GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label)
311 fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)nClusters);
314 } else { // by comparing cluster/TrackRef positions
316 // look for the corresponding simulated track if any
317 TIter next(trackRefStore.CreateIterator());
318 AliMUONTrack* trackRef;
319 while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
321 // check compatibility
323 if (TrackMatched(muonTrack, *trackRef, f, sigmaCut)) {
324 matchedTrackRef = trackRef;
325 fractionOfMatchCluster = f;
333 return matchedTrackRef;