Fixes for #86059: Install data when ALICE_ROOT does not point to source (Christian)
[u/mrichter/AliRoot.git] / MUON / DIMUONFakes.C
CommitLineData
30b3fe0f 1#if !defined(__CINT__) || defined(__MAKECINT__)
2// ROOT includes
30b3fe0f 3#include <TFile.h>
4#include <TH1.h>
5#include <TCanvas.h>
6#include <Riostream.h>
7#include <TROOT.h>
8#include <TClonesArray.h>
30b3fe0f 9#include <TLorentzVector.h>
10
11// STEER includes
12#include "AliLog.h"
30b3fe0f 13#include "AliESDEvent.h"
14#include "AliESDMuonTrack.h"
30b3fe0f 15#include "AliCDBManager.h"
16
17// MUON includes
a99c3449 18#include "AliMUONCDB.h"
30b3fe0f 19#include "AliMUONTrack.h"
20#include "AliMUONVTrackStore.h"
21#include "AliMUONTrackParam.h"
30b3fe0f 22#include "AliMUONRecoCheck.h"
23#include "AliMUONVCluster.h"
24#include "AliMUONRecoParam.h"
25#endif
26
27/// \ingroup macros
28/// \file DIMUONFakes.C
29///
30/// \author Ph. Pillot, Subatech, March. 2009
31///
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
35
f202486b 36Double_t sigmaCut = -1.;
30b3fe0f 37
38//-----------------------------------------------------------------------
39void DIMUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
a99c3449 40 const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/",
41 const TString ocdbPath = "local://$ALICE_ROOT/OCDB")
30b3fe0f 42{
43
44 //Reset ROOT and connect tree file
45 gROOT->Reset();
46
47 // File for histograms and histogram booking
48 TFile *histoFile = new TFile("DiFakes.root", "RECREATE");
49
23035434 50 TH1F *hMass = new TH1F("hMass", "Dimuon mass distribution (GeV/c^{2})", 100, 0., 12.);
30b3fe0f 51 TH1F *hMassM = new TH1F("hMassM", "matched track mass distribution (GeV/c^{2})", 100, 0., 12.);
52 TH1F *hMassF = new TH1F("hMassF", "fake track mass distribution (GeV/c^{2})", 100, 0., 12.);
23035434 53 TH1F *hP = new TH1F("hP", "Dimuon P distribution (GeV/c)", 100, 0., 200.);
30b3fe0f 54 TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
55 TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
23035434 56 TH1F *hPt = new TH1F("hPt", "Dimuon Pt distribution (GeV/c)", 100, 0., 20.);
30b3fe0f 57 TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
58 TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
23035434 59 TH1F *hY = new TH1F("hY"," Dimuon rapidity distribution",100,-10,0);
30b3fe0f 60 TH1F *hYM = new TH1F("hYM"," matched track rapidity distribution",100,-10,0);
61 TH1F *hYF = new TH1F("hYF"," fake track rapidity distribution",100,-10,0);
23035434 62 TH1F *hEta = new TH1F("hEta"," Dimuon pseudo-rapidity distribution",100,-10,0);
30b3fe0f 63 TH1F *hEtaM = new TH1F("hEtaM"," matched track pseudo-rapidity distribution",100,-10,0);
64 TH1F *hEtaF = new TH1F("hEtaF"," fake track pseudo-rapidity distribution",100,-10,0);
23035434 65 TH1F *hPhi = new TH1F("hPhi"," Dimuon phi distribution",100,-1.,9.);
30b3fe0f 66 TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.);
67 TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.);
68
a99c3449 69 // link to reconstructed and simulated tracks
70 AliMUONRecoCheck rc(esdFileName, SimDir);
30b3fe0f 71
a99c3449 72 // load necessary data from OCDB
73 AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
74 AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
75 AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
76 if (!recoParam) return;
30b3fe0f 77
a99c3449 78 // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
f202486b 79 sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
30b3fe0f 80
81 TLorentzVector vMu1, vMu2, vDiMu;
82
83 // Loop over ESD events
84 FirstEvent = TMath::Max(0, FirstEvent);
a99c3449 85 LastEvent = (LastEvent>=0) ? TMath::Min(rc.NumberOfEvents() - 1, LastEvent) : rc.NumberOfEvents() - 1;
30b3fe0f 86 for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
87
a99c3449 88 // get reconstructed and simulated tracks
89 AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(iEvent, kFALSE);
30b3fe0f 90 AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
a99c3449 91 if (!muonTrackStore || !trackRefStore) continue;
30b3fe0f 92
93 // loop over ESD tracks and flag them
a99c3449 94 const AliESDEvent* esd = rc.GetESDEvent();
30b3fe0f 95 Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
96 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
97
a99c3449 98 AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
30b3fe0f 99
100 // skip ghosts
a99c3449 101 if (!esdTrack->ContainTrackerData()) continue;
102
103 // find the corresponding MUON track
104 AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID());
30b3fe0f 105
106 // try to match the reconstructed track with a simulated one
f202486b 107 Int_t nMatchClusters = 0;
108 AliMUONTrack* matchedTrackRef = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClusters, useLabel, sigmaCut);
30b3fe0f 109
110 // take actions according to matching result
111 if (matchedTrackRef) {
112
113 // flag matched tracks
a99c3449 114 esdTrack->SetLabel(matchedTrackRef->GetUniqueID());
30b3fe0f 115
116 // remove already matched trackRefs
117 trackRefStore->Remove(*matchedTrackRef);
118
119 } else {
120
121 // flag fake tracks
a99c3449 122 esdTrack->SetLabel(-1);
30b3fe0f 123
124 }
125
126 }
127
128 // double loop over ESD tracks, build pairs and fill histograms according to their label
129 for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
130 AliESDMuonTrack* muonTrack1 = esd->GetMuonTrack(iTrack1);
131
132 // skip ghosts
133 if (!muonTrack1->ContainTrackerData()) continue;
134
135 // get track info
136 Short_t charge1 = muonTrack1->Charge();
137 Int_t label1 = muonTrack1->GetLabel();
138 muonTrack1->LorentzP(vMu1);
139
140 for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
141 AliESDMuonTrack* muonTrack2 = esd->GetMuonTrack(iTrack2);
142
143 // skip ghosts
144 if (!muonTrack2->ContainTrackerData()) continue;
145
146 // keep only opposite sign pairs
147 Short_t charge2 = muonTrack2->Charge();
148 if (charge1*charge2 > 0) continue;
149
150 // get track info
151 Int_t label2 = muonTrack2->GetLabel();
152 muonTrack2->LorentzP(vMu2);
153
154 // compute kinematics of the pair
155 vDiMu = vMu1 + vMu2;
156 Float_t mass = vDiMu.M();
157 Float_t p = vDiMu.P();
158 Float_t pt = vDiMu.Pt();
159 Float_t y = vDiMu.Rapidity();
160 Float_t eta = vDiMu.Eta();
23035434 161 Float_t phi = vDiMu.Phi();
30b3fe0f 162 if (phi < 0) phi += 2.*TMath::Pi();
163
164 // fill global histograms
165 hMass->Fill(mass);
166 hP->Fill(p);
167 hPt->Fill(pt);
168 hY->Fill(y);
169 hEta->Fill(eta);
170 hPhi->Fill(phi);
171
172 // fill histograms according to labels
173 if (label1 >= 0 && label2 >= 0) {
174
175 hMassM->Fill(mass);
176 hPM->Fill(p);
177 hPtM->Fill(pt);
178 hYM->Fill(y);
179 hEtaM->Fill(eta);
180 hPhiM->Fill(phi);
181
182 } else {
183
184 hMassF->Fill(mass);
185 hPF->Fill(p);
186 hPtF->Fill(pt);
187 hYF->Fill(y);
188 hEtaF->Fill(eta);
189 hPhiF->Fill(phi);
190
191 }
192
193 }
194
195 }
196
197 } // end of loop over events
198
199 // plot results
200 TCanvas cDiFakesSummary("cDiFakesSummary","cDiFakesSummary",900,600);
201 cDiFakesSummary.Divide(3,2);
202 cDiFakesSummary.cd(1);
203 cDiFakesSummary.GetPad(1)->SetLogy();
204 hMass->Draw();
205 hMass->SetMinimum(0.5);
206 hMassM->Draw("same");
207 hMassM->SetLineColor(4);
208 hMassF->Draw("same");
209 hMassF->SetLineColor(2);
210 hMassF->SetFillColor(2);
211 hMassF->SetFillStyle(3017);
212 cDiFakesSummary.cd(2);
213 cDiFakesSummary.GetPad(3)->SetLogy();
214 hP->Draw();
215 hP->SetMinimum(0.5);
216 hPM->Draw("same");
217 hPM->SetLineColor(4);
218 hPF->Draw("same");
219 hPF->SetLineColor(2);
220 hPF->SetFillColor(2);
221 hPF->SetFillStyle(3017);
222 cDiFakesSummary.cd(3);
223 cDiFakesSummary.GetPad(4)->SetLogy();
224 hPt->Draw();
225 hPt->SetMinimum(0.5);
226 hPtM->Draw("same");
227 hPtM->SetLineColor(4);
228 hPtF->Draw("same");
229 hPtF->SetLineColor(2);
230 hPtF->SetFillColor(2);
231 hPtF->SetFillStyle(3017);
232 cDiFakesSummary.cd(4);
233 cDiFakesSummary.GetPad(2)->SetLogy();
234 hY->Draw();
235 hY->SetMinimum(0.5);
236 hYM->Draw("same");
237 hYM->SetLineColor(4);
238 hYF->Draw("same");
239 hYF->SetLineColor(2);
240 hYF->SetFillColor(2);
241 hYF->SetFillStyle(3017);
242 cDiFakesSummary.cd(5);
243 cDiFakesSummary.GetPad(5)->SetLogy();
244 hEta->Draw();
245 hEta->SetMinimum(0.5);
246 hEtaM->Draw("same");
247 hEtaM->SetLineColor(4);
248 hEtaF->Draw("same");
249 hEtaF->SetLineColor(2);
250 hEtaF->SetFillColor(2);
251 hEtaF->SetFillStyle(3017);
252 cDiFakesSummary.cd(6);
253 cDiFakesSummary.GetPad(6)->SetLogy();
254 hPhi->Draw();
255 hPhi->SetMinimum(0.5);
256 hPhiM->Draw("same");
257 hPhiM->SetLineColor(4);
258 hPhiF->Draw("same");
259 hPhiF->SetLineColor(2);
260 hPhiF->SetFillColor(2);
261 hPhiF->SetFillStyle(3017);
262
263 // save results
264 histoFile->cd();
265 histoFile->Write();
266 cDiFakesSummary.Write();
267 histoFile->Close();
268
30b3fe0f 269}
270