]>
Commit | Line | Data |
---|---|---|
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 | 36 | Double_t sigmaCut = -1.; |
30b3fe0f | 37 | |
38 | //----------------------------------------------------------------------- | |
39 | void 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 | ||
590a335d | 47 | AliLog::SetClassDebugLevel("AliMCEvent",-1); |
48 | ||
30b3fe0f | 49 | // File for histograms and histogram booking |
50 | TFile *histoFile = new TFile("DiFakes.root", "RECREATE"); | |
51 | ||
23035434 | 52 | TH1F *hMass = new TH1F("hMass", "Dimuon mass distribution (GeV/c^{2})", 100, 0., 12.); |
30b3fe0f | 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.); | |
23035434 | 55 | TH1F *hP = new TH1F("hP", "Dimuon P distribution (GeV/c)", 100, 0., 200.); |
30b3fe0f | 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.); | |
23035434 | 58 | TH1F *hPt = new TH1F("hPt", "Dimuon Pt distribution (GeV/c)", 100, 0., 20.); |
30b3fe0f | 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.); | |
23035434 | 61 | TH1F *hY = new TH1F("hY"," Dimuon rapidity distribution",100,-10,0); |
30b3fe0f | 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); | |
23035434 | 64 | TH1F *hEta = new TH1F("hEta"," Dimuon pseudo-rapidity distribution",100,-10,0); |
30b3fe0f | 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); | |
23035434 | 67 | TH1F *hPhi = new TH1F("hPhi"," Dimuon phi distribution",100,-1.,9.); |
30b3fe0f | 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.); | |
70 | ||
a99c3449 | 71 | // link to reconstructed and simulated tracks |
72 | AliMUONRecoCheck rc(esdFileName, SimDir); | |
30b3fe0f | 73 | |
a99c3449 | 74 | // load necessary data from OCDB |
75 | AliCDBManager::Instance()->SetDefaultStorage(ocdbPath); | |
590a335d | 76 | AliCDBManager::Instance()->SetSpecificStorage("GRP/GRP/Data","local://."); |
a99c3449 | 77 | AliCDBManager::Instance()->SetRun(rc.GetRunNumber()); |
78 | AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam(); | |
79 | if (!recoParam) return; | |
30b3fe0f | 80 | |
a99c3449 | 81 | // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used |
f202486b | 82 | sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking(); |
590a335d | 83 | // compute the mask of requested stations from recoParam |
84 | UInt_t requestedStationMask = 0; | |
85 | for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) requestedStationMask |= ( 1 << i ); | |
86 | // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible | |
87 | Bool_t request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates(); | |
30b3fe0f | 88 | |
89 | TLorentzVector vMu1, vMu2, vDiMu; | |
590a335d | 90 | Int_t nReconstructiblePairs = 0; |
91 | Int_t nReconstructedPairs = 0; | |
92 | Int_t nMatchedPairs = 0; | |
30b3fe0f | 93 | |
94 | // Loop over ESD events | |
95 | FirstEvent = TMath::Max(0, FirstEvent); | |
a99c3449 | 96 | LastEvent = (LastEvent>=0) ? TMath::Min(rc.NumberOfEvents() - 1, LastEvent) : rc.NumberOfEvents() - 1; |
30b3fe0f | 97 | for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) { |
98 | ||
a99c3449 | 99 | // get reconstructed and simulated tracks |
100 | AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(iEvent, kFALSE); | |
30b3fe0f | 101 | AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent); |
a99c3449 | 102 | if (!muonTrackStore || !trackRefStore) continue; |
30b3fe0f | 103 | |
590a335d | 104 | // count the number of reconstructible pairs |
105 | Int_t nMuPlus = 0, nMuMinus = 0; | |
106 | TIter next(trackRefStore->CreateIterator()); | |
107 | AliMUONTrack* trackRef; | |
108 | while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) { | |
109 | if (trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) { | |
110 | if (trackRef->GetTrackParamAtVertex()->GetCharge() > 0) nMuPlus++; | |
111 | else nMuMinus++; | |
112 | } | |
113 | } | |
114 | nReconstructiblePairs += nMuPlus*nMuMinus; | |
115 | ||
30b3fe0f | 116 | // loop over ESD tracks and flag them |
a99c3449 | 117 | const AliESDEvent* esd = rc.GetESDEvent(); |
30b3fe0f | 118 | Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ; |
119 | for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { | |
120 | ||
a99c3449 | 121 | AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack); |
30b3fe0f | 122 | |
123 | // skip ghosts | |
a99c3449 | 124 | if (!esdTrack->ContainTrackerData()) continue; |
125 | ||
126 | // find the corresponding MUON track | |
127 | AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID()); | |
30b3fe0f | 128 | |
129 | // try to match the reconstructed track with a simulated one | |
f202486b | 130 | Int_t nMatchClusters = 0; |
131 | AliMUONTrack* matchedTrackRef = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClusters, useLabel, sigmaCut); | |
30b3fe0f | 132 | |
133 | // take actions according to matching result | |
134 | if (matchedTrackRef) { | |
135 | ||
136 | // flag matched tracks | |
a99c3449 | 137 | esdTrack->SetLabel(matchedTrackRef->GetUniqueID()); |
30b3fe0f | 138 | |
139 | // remove already matched trackRefs | |
140 | trackRefStore->Remove(*matchedTrackRef); | |
141 | ||
142 | } else { | |
143 | ||
144 | // flag fake tracks | |
a99c3449 | 145 | esdTrack->SetLabel(-1); |
30b3fe0f | 146 | |
147 | } | |
148 | ||
149 | } | |
150 | ||
151 | // double loop over ESD tracks, build pairs and fill histograms according to their label | |
152 | for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) { | |
153 | AliESDMuonTrack* muonTrack1 = esd->GetMuonTrack(iTrack1); | |
154 | ||
155 | // skip ghosts | |
156 | if (!muonTrack1->ContainTrackerData()) continue; | |
157 | ||
158 | // get track info | |
159 | Short_t charge1 = muonTrack1->Charge(); | |
160 | Int_t label1 = muonTrack1->GetLabel(); | |
161 | muonTrack1->LorentzP(vMu1); | |
162 | ||
163 | for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) { | |
164 | AliESDMuonTrack* muonTrack2 = esd->GetMuonTrack(iTrack2); | |
165 | ||
166 | // skip ghosts | |
167 | if (!muonTrack2->ContainTrackerData()) continue; | |
168 | ||
169 | // keep only opposite sign pairs | |
170 | Short_t charge2 = muonTrack2->Charge(); | |
171 | if (charge1*charge2 > 0) continue; | |
172 | ||
590a335d | 173 | nReconstructedPairs++; |
174 | ||
30b3fe0f | 175 | // get track info |
176 | Int_t label2 = muonTrack2->GetLabel(); | |
177 | muonTrack2->LorentzP(vMu2); | |
178 | ||
179 | // compute kinematics of the pair | |
180 | vDiMu = vMu1 + vMu2; | |
181 | Float_t mass = vDiMu.M(); | |
182 | Float_t p = vDiMu.P(); | |
183 | Float_t pt = vDiMu.Pt(); | |
184 | Float_t y = vDiMu.Rapidity(); | |
185 | Float_t eta = vDiMu.Eta(); | |
23035434 | 186 | Float_t phi = vDiMu.Phi(); |
30b3fe0f | 187 | if (phi < 0) phi += 2.*TMath::Pi(); |
188 | ||
189 | // fill global histograms | |
190 | hMass->Fill(mass); | |
191 | hP->Fill(p); | |
192 | hPt->Fill(pt); | |
193 | hY->Fill(y); | |
194 | hEta->Fill(eta); | |
195 | hPhi->Fill(phi); | |
196 | ||
197 | // fill histograms according to labels | |
198 | if (label1 >= 0 && label2 >= 0) { | |
199 | ||
590a335d | 200 | nMatchedPairs++; |
201 | ||
30b3fe0f | 202 | hMassM->Fill(mass); |
203 | hPM->Fill(p); | |
204 | hPtM->Fill(pt); | |
205 | hYM->Fill(y); | |
206 | hEtaM->Fill(eta); | |
207 | hPhiM->Fill(phi); | |
208 | ||
209 | } else { | |
210 | ||
211 | hMassF->Fill(mass); | |
212 | hPF->Fill(p); | |
213 | hPtF->Fill(pt); | |
214 | hYF->Fill(y); | |
215 | hEtaF->Fill(eta); | |
216 | hPhiF->Fill(phi); | |
217 | ||
218 | } | |
219 | ||
220 | } | |
221 | ||
222 | } | |
223 | ||
224 | } // end of loop over events | |
225 | ||
226 | // plot results | |
227 | TCanvas cDiFakesSummary("cDiFakesSummary","cDiFakesSummary",900,600); | |
228 | cDiFakesSummary.Divide(3,2); | |
229 | cDiFakesSummary.cd(1); | |
230 | cDiFakesSummary.GetPad(1)->SetLogy(); | |
231 | hMass->Draw(); | |
232 | hMass->SetMinimum(0.5); | |
233 | hMassM->Draw("same"); | |
234 | hMassM->SetLineColor(4); | |
235 | hMassF->Draw("same"); | |
236 | hMassF->SetLineColor(2); | |
237 | hMassF->SetFillColor(2); | |
238 | hMassF->SetFillStyle(3017); | |
239 | cDiFakesSummary.cd(2); | |
240 | cDiFakesSummary.GetPad(3)->SetLogy(); | |
241 | hP->Draw(); | |
242 | hP->SetMinimum(0.5); | |
243 | hPM->Draw("same"); | |
244 | hPM->SetLineColor(4); | |
245 | hPF->Draw("same"); | |
246 | hPF->SetLineColor(2); | |
247 | hPF->SetFillColor(2); | |
248 | hPF->SetFillStyle(3017); | |
249 | cDiFakesSummary.cd(3); | |
250 | cDiFakesSummary.GetPad(4)->SetLogy(); | |
251 | hPt->Draw(); | |
252 | hPt->SetMinimum(0.5); | |
253 | hPtM->Draw("same"); | |
254 | hPtM->SetLineColor(4); | |
255 | hPtF->Draw("same"); | |
256 | hPtF->SetLineColor(2); | |
257 | hPtF->SetFillColor(2); | |
258 | hPtF->SetFillStyle(3017); | |
259 | cDiFakesSummary.cd(4); | |
260 | cDiFakesSummary.GetPad(2)->SetLogy(); | |
261 | hY->Draw(); | |
262 | hY->SetMinimum(0.5); | |
263 | hYM->Draw("same"); | |
264 | hYM->SetLineColor(4); | |
265 | hYF->Draw("same"); | |
266 | hYF->SetLineColor(2); | |
267 | hYF->SetFillColor(2); | |
268 | hYF->SetFillStyle(3017); | |
269 | cDiFakesSummary.cd(5); | |
270 | cDiFakesSummary.GetPad(5)->SetLogy(); | |
271 | hEta->Draw(); | |
272 | hEta->SetMinimum(0.5); | |
273 | hEtaM->Draw("same"); | |
274 | hEtaM->SetLineColor(4); | |
275 | hEtaF->Draw("same"); | |
276 | hEtaF->SetLineColor(2); | |
277 | hEtaF->SetFillColor(2); | |
278 | hEtaF->SetFillStyle(3017); | |
279 | cDiFakesSummary.cd(6); | |
280 | cDiFakesSummary.GetPad(6)->SetLogy(); | |
281 | hPhi->Draw(); | |
282 | hPhi->SetMinimum(0.5); | |
283 | hPhiM->Draw("same"); | |
284 | hPhiM->SetLineColor(4); | |
285 | hPhiF->Draw("same"); | |
286 | hPhiF->SetLineColor(2); | |
287 | hPhiF->SetFillColor(2); | |
288 | hPhiF->SetFillStyle(3017); | |
289 | ||
290 | // save results | |
291 | histoFile->cd(); | |
292 | histoFile->Write(); | |
293 | cDiFakesSummary.Write(); | |
294 | histoFile->Close(); | |
295 | ||
590a335d | 296 | // print results |
297 | printf("\n"); | |
298 | printf("nb of reconstructible OS pairs: %d \n", nReconstructiblePairs); | |
299 | printf("nb of reconstructed OS pairs: %d \n", nReconstructedPairs); | |
300 | printf("nb of reconstructed OS pairs matched with trackRefs: %d \n", nMatchedPairs); | |
301 | printf("\n"); | |
302 | ||
30b3fe0f | 303 | } |
304 |